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ABSTRACT 


One experimental technique based on the Faraday rotation effect of 
radio waves is presented for measuring electron density in the nighttime 
ionosphere at midlatitudes. In the experiment high frequency linearly- 
polarized radio signals are transmitted to a linearly-polarized receiving 
system located in a spinning rocket moving through the ionosphere. 

Faraday rotation is observed in the reference plane of the rocket as a 
change in frequency of the detected receiver output. Computer analysis 
measures this frequency change and uses the information to obtain electron- 
density data. System performance is evaluated and some sources of error 
are identified. It is found that the dominant controlling factor for 
accuracy is signal-to-noise ratio. The data obtained from the above method 
has been useful in calibrating a Langmuir probe experiment for electron- 
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density values of 100 cm and greater. Data from two rocket flights are 
presented to illustrate the experiment. 
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1. INTRODUCTION 

The University of Illinois rocket-radio propagation experiment was 
designed in 1963 and first used for the International Quiet Sun Year (IQSY) 
in 1964 and 1965. The propagation experiment is used in conjunction with a 
Langmuir probe to provide a detailed description of electron concentration 
versus altitude in the lower ionosphere. The Langmuir probe is able to follow 
the detailed variation (fine structure) of the ionosphere's electron concen- 
tration profile, but is unable to convert its measurement directly into 
electron density of sufficient accuracy. Therefore, the propagation experi- 
ment is required to provide an independent measure of average electron 
concentration that can be used tc calibrate the more detailed measurements of 
the Langmuir probe. 

In general, it can be stated that for daytime launches at midlatitudes 
the propagation experiment, used since 1964 with only slight modification, 
has been very successful in providing the necessary calibration. For night- 
time launches at midlatitudes, however, the propagation experiment has not 
been satisfactory. To overcome this problem and to help identify possible 
limiting sources of error, the original experiment has been modified and 
simplified. The original experiment measured both differential absorption 
and Faraday rotation; the modified experiment measures Faraday rotation only. 

The modified experiment was first used for the flight of Nike Apache 
14.534, launched on 10 August 1977. It was subsequently used for the flight 
of Nike Apache 14.533, launched on 5 January 1978. Both flights used the 
modified experiment at a frequency at 2.225 MHz. The flights also used the 
original experiment at 3.385 MHz in order to allow a comparison of the two 
experiments. The data from these flights are presented in this report. 

The organization of this report will be to first present, in Chapter 2, 
the reasons for making measurements of electron density in the nighttime as 
well as a comparison between day and night profiles. The principle of 
Faraday rotation measurement and its application to both the original and 
modified experiment are presented in Chapter 3. Chapter 4 contains a 
description of the modified experiment system including the ground-based 
transmitter, the rocket payload and the necessary telemetry and data 
recording support. The data reduction process used in the experiment is 
presented in Chapter 5 where the final electron-density profiles are shown. 

The experimental limitations that have been identified in both the hardware 
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and software are discussed in Chapter 6. Finally, Chapter 7 contains a 
summary and makes conclusions about the modified propagation ex^riment . In 
addition, recommendations for future work are given. 

It will be seen that both the original and modified experiments are 
capable of measuring electron densities as small as 100 cm” J . The modified 
experiment is, however, easier to implement and is preferred for nighttime 
measurements . 
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2. THE NIGHTTIME # REGION 

2.1 General Properties cold Reasons for Measuring Nighttime Profiles 

The ionosphere plays an important role in a variety of effects in radio 
propagation [Laurence et al. , 1964] and in upper atmospheric physics 
[Ratoliffe, 1972]. Electron density and its altitude distribution are 
essential information in the analysis of these effects. As a ^ . 

electron-density profile measurements have been of fundamenta .mporta..' 
Further, considerable effort has been made in identifying and studying the 
sources of ionization production and the role of atmospheric dynamic? as 
they affect these profiles. This effort has been quite successful in 
explaining the major features of daytime profiles, but work for nighttime 
profiles still remains. 

The major hindrance has been the difficulty of obtaining detailed and 
accurate nighttime electron-density profiles. In addition to their importance 
in the analysis of other physical effects, these profiles can be used in 
revealing information about the ionosphere itself. Smith et al. [1974] used 
nighttime profiles to calculate ionization rates in the ionosphere, while 
Vo 8 8 and Smith [1977] correlated certain profile characteristics with geo- 
magnetic activity to identify energetic particles as a major source of 
nighttime ionization in the 120 to 200 km region under disturbed geomagnetic 
conditions. 

A daytime ionosphere is shown in Figure 2.1 (profile 1). Electron 

densities in the daytime are determined by the solar radiation in the UV and 

X-ray parts of the spectrum. Typical electron densities are on the order of 
5 -3 

10 cm and are reproc jcible from day to day for a given solar zenith angle. 
One exception to this is the sporadic-# layers which can occur in both day 
and night profiles in the 90 to 110 km region. Figure 2.1 shows sporadic-# 
layers in profiles 1 and 2 at 100 and 106 km, respectively. 

After sunset electron density falls off rapidly and within two hours 

3 -3 

reaches the typical nighttime value of 10 cm . The nighttime ionosphere 
has been shown to have substantial amounts of temporal variations and complex 
altitude structure [Smith, 1970]. In addition, very little radiowave 
absorption is observed in the night [Diokinson et al., 1976]. Figure 2.1 
(profile 2) shows a nighttime electron-density profile taken near local 
midnight which demonstrates the more compiex structure. As can be seen, in 
addition to the intense sporadic-# layer at 106 km a much wider layer is 



ALTITUDE (km) 


©NIKE APACHE 14.514 
10 AUGUST, 1973 
1005 EST 
WALLOPS ISLANO 

(D NIKE APACHE 14.522 
29 JUNE , 1974 
2340 EST 
WALLOPS ISLANO 


ELECTRON DENSITY (cm' 3 ) 


Figure 2.1 Profile 1 is representative of the daytime electron 
density (10 3 cm" 3 ) . Towards night this electron 
density rapidly decays to a lower nighttime equilib- 
rium level of 10 3 cm' 3 (profile 2) . The nighttime 
density profile shows much variability and structure 
[Voss and Smith, 1977] . 
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present between 110 km and 17S km; this is the intermediate layer. The 
intermediate layer is located between the E and F layers with a vertical 
extent, typically, of tens of kilometers and is the result of a wind shear. 
Through the night this layer tends to move down in altitude following the 
atmospheric tide. The valleys below and above the layer may have electron 

_3 

concentrations as low as 200 cm . 

2.2 Methods of Nighttime Profile Measurements using Rockets 

Knight [1972] has nude a survey of published nighttime electron-density 
profiles. Of the profiles that were obtained using rocket measurements, most 
use versions of the Langmuir ; ;be. However, the techniques used to calibrate 
the probe varied. By far the simplest procedure is to rely solely on 
ground-based experiments such as ionosondes. A second class of experiments 
is that which used an independent probing system to measure the electron 
density; for example, the resonance probe [Mirao and Muraoka , 1967] which 
measures plasma frequency, and the standing wave impedance probe [Ultfick 
et al.y 1967] which measures the input impedance of an antenna deployed in 
the ionospheric plasma. 

A third class of experiments is that using rocket-borne radio propaga- 
tion experiments. The major advantage of •'•his method is that the radio 
signal propagates through a region of the ionosphere that has been unaffected 
by the passage of the rocket except for a localized region around the rocket 
[. Bowhilly 1964]. Within this classification a variety of techniques have 
been used. In general, the radio propagatior experiments can be divided into 
two categories according to which property of the propagated wave is measured: 
absorption and phase (or frequency!. Many experimental systems incorporate 
methods to make measurements in both categories. The absorption measurement 
provides information on the imaginary part of the complex index of refraction 
for the ionosphere. Similarly, phase or frequency measurement yields informa- 
tion on the real part of the index of refraction. Both the real and 
imaginary parts can be used to calculate electron density. Further subdivi- 
sion within each category is possible based upon the frequency of the 
propagated wave: VLF experiments have been reported that operate below the 

ionospheric, plasma and gyrofroquencies; there are also HF and VHF experiments 
operating abuve these frequencies. Examples ot these categories applied to 
the nighttime ionosphere will now be given. 
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Cartwright [1964] and Kimura et al. [1967] have demonstrated the use 
of VLF and LF Doppler frequency measurement techniques at 10.2 and 40 kHz, 
respectively. Each used a remotely situated transmitter with the signal 
propagating in the earth-ionosphere waveguide. By comparing the Doppler- 
shifted received signal with a stable local oscillator on board the rocket 
and telemetering the beat frequency to the ground, the Doppler shift is 
measured. Because Doppler shift is affected by the index of refraction of 
the medium, electron density can then be extracted. To amplify the Doppler 
shift before measurement, Kimura first multiplied the beat frequency by 200 
thus increasing the effect. Cartwright's final electron-density profile is 
shown in Figure 2.2. This profile is significant as it is the first in situ 
observation of the intermediate layer in a nighttime electron-density profile. 

A dispersive Doppler technique at much greater frequencies was used by 
Aikin and Blumle [1968] to measure nighttime profiles at the equator. Their 
experiment was based on the technique first introduced by Seddon [1953J. The 
experiment uses two harmonically related frequencies (at 24.5 and 73.5 MHz) 
transmitted from the rocket. These signals are received on the ground where 
the low frequency is multiplied up to the higher and their beat frequency 
determined. This frequency contains index of refraction information which 
can be used to give electron density. Bowhill [1964] gives a detailed 
description of this experimental technique. 

Nighttime techniques to measure absorption have in general not been 
successful in obtaining accurate electron-density profiles. This is because 
of the low level of absorption present in the nighttime ionosphere and is 
further complicated by the presence of ionospheric reflections causing 
standing waves. Some success was reported by Hall and Books [1965] using 
a ground-based LF transmitter at 202 kHz and a rocket payload designed to 
sense the field strength of the signal. However, the presence of standing 
waves reduced the experiment's accuracy. 

2.3 Conclusion 

This chapter has discussed the characteristics of the nighttime iono- 
sphere and presented some possible methods to measure the electron-density 
profiles. Experience with the rocket program of the Aeronomy Laboratory has 
shown that a Langmuir probe combined with a radio propagation experiment is 
very effective. As has been mentioned earlier, electron-density measurement 
using radio propagation is in general more difficult at night than during the 
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Figure 2.2 Nighttime electron-density profile as 

measured by Cartwright [1964] using VLF 
Doppler experiment. An intermediate 
layer is shown at 140 km. 
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day. The reasons for this can now be summarized. Firstly, the nightime 
electron density is typically two orders of magnitude smaller than the 
daytime value, which causes a corresponding reduction in all radio iono- 
spheric effects. Secondly, the more complex structure and regions of low 
electron density which occur in the nighttime E region require that the 
experiment has good altitude resolution and sensitivity. Finally, because 
of the possibility of strong reflections (in the virtual absence of 
absorption) any experiment must either be insensitive to them or account 
for their existence in the analysis of the results. The following chapters 
will present a relatively simple radio propagation experiment which appears 
to satisfy these nighttime requirements. 
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3. THE FARADAY ROTATION EXPERIMENT 

3.1 Introduction 

The original propagation experiment and the modified version for night- 
time use, to be described in detail in a later chapter, both take advantage 
of the fact that a radio wave propagating through the ionosphere is modified 
in characteristic ways. If a rocket carrying a receiver can detect how 
the radio signal has been altered, information on electron concentration 
can be determined as a function of altitude. This chapter will discuss 
the use of Faraday rotation as a technique for measuring electron 
concentration. 

Faraday rotation is a phenomenon in which a linearly polarized signal 
has its direction of polarization rotated as it travels. The use of this 
effect as a measurement technique for the ionosphere can be considered to 
be part of a general class of experiments designed to measure differential 
phase velocities between differing plane wave propagating modes or 
signals of different frequencies [ Ratcliffe , 1972], 

The chapter begins with a section discussing the general theory of 
Faraday rotation. This is followed by two sections describing how Faraday 
rotation is measured in both the original and the modified radio propagation 
experiments. Additionally, the section dealing with the original experiment 
contains a brief system description to provide a basis for comparing it with 
the modified experiment. 

Although the particulars of the modified experiment are to be discussed 
in the following chapter, it will be useful at this point to briefly discuss 
the reasoning involved in modifying the original experiment, as this will 
illuminate why the particular choice of transmitted radio signal used in the 
modified experiment was made. As has been previously mentioned, the original 
propagation experiment had not been successful in providing usable nighttime 
data; this was because the errors internal to the experimental system were 
larger than the effects to be measured. However, because of the system's 
complexity, any significant sources of error could not be positively 
identified. It was felt that if the system's complexity could be reduced, 
some sources of error might be eliminated while simplifying the analysis of 
any remaining error sources. Further, because of time constraints prior to 
the launch of Nike Apache 14.534, no substantial hardware development for a 
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new experiment could be completed in time for this rocket flight. Therefore, 
it was concluded that the design of the new experiment must be based on 
existing equipment. Also, because the original experiment would still be 
used for daytime flights, this equipment could not be extensively modified. 

There are two general areas that introduce operational complexity into 
the original experiment: these are the feedback amplitude control loop and 

its requirement for independent generation of both propagating modes. These 
are used in the measurement of differential absorption, which in the night* 
time is negligible. Therefore, by not attempting to measure differential 
absorption and relying solely on Faraday rotation information, the related 
hardware for these two areas could be eliminated, greatly simplifying the 
system’s operation. It will be shown that this new experiment can success- 
fully measure the electron-density profile of the nighttime E region. 

3.2 Faraday Rotation 

Michael Faraday first demonstrated (in 1845) that the plane of polari- 
zation of a light beam is rotated when passed through certain substances 
having a magnetic field aligned along the light's direction of propagation. 

It was further found that the amount of rotation was proportional to the 
magnitude of the magnetic field and the distance traversed in the substance. 
This phenomenon of polarization rotation has since acquired the name of 
Faraday rotation or Faraday effect. 

For the purposes of this report we are interested in describing the 
.araday rotation of a radio wave in the lower ionosphere. To do this the 
necessary theoretical development begins with the work of Appleton [1932] 
and t.artree [1931], By using the results of their work ionospheric Faraday 
rotation and the role electron concentration plays in this effect can be 
described. 

The presence of a magnetic field in a region of free electrons, such 
the ionosphere, causes a plane electromagnetic wave to split into two 
characteristic modes having different polarizations and indices of refraction. 
The Applet on -Hart ree theory can be used to determine both the polarization and 
index of refraction for each mode. Their theory is based upon the following 
assumptions: the ionospheric medium is electrically neutral, electron- 

neutral collision frequency is independent of electron energy, the 
magnetic field is uniform, and a background of essentially stationary 
ions exists. 
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The Appleton-Hartree analysis begins by assuming that a plane electro- 
magnetic wave (satisfying Maxwell's equations) is propagating in a positive 
z direction 

h,.t (3.1) 


where k is the propagation constant and u is the frequency of the radio 
signal. I is a complex vector containing amplitude and phase information 
for the directional components of the electric field 

^ A A 

E = x E + y E (3.2) 

X y 

Figure 3.1 presents a simplified diagram depicting the situation. Also 
indicated is the earth's magnetic field, B, with an assumed direction, and 
the ionosphere with electron density, N . 

The motion of ionospheric free electrons is affected both by the earth's 
magnetic fiijld and by the travelling electromagnetic wave. This is reflected 
in the equation of motion for each electron 
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where r is the electron displacement from some rest position, v is the *!cctron- 
neutral collision frequency and m and e are the mass and charge of the electron. 
The above equation makes the assumption that the magnetic field associated 
with the travelling wave is very much smaller than the earth's magnetic field. 

By requiring that an electromagnetic wave simultaneously satisfy all 
Maxwell's equations as well as equation (3.3) the work of Appleton and of 
Hartree demonstrated that only two possible travelling waves of the form in 
equation (3.1) could propagate through the ionosphere. Each of these has a 
unique polarization and index of refraction; before these can be given the 
following standard symbols must be introduced: 


X 



4. 

U) 


N e 
e 


me w 
o 


(3.4) 


= dJL 


CO 


mto 





(3.5) 



12 



figure 3.1 Simplified model at midlatitudes (Wallops Island) for plane- 
wave propagation through a region with electron density, PJ e , 
and magnetic field, S. At midlatitudes the earth’s magnetic 
field is nearly antiparallel to the z axis. 
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Z 



(3.6) 


where is the plasma frequency and u> ff is the gyrofrequency. Also, 

Tj. = 7 cose and Y T - Y sine, where 8 is the angle between the magnetic 
field and the direction of propagation. 

Using the above defined terms, the allowed propagating modes have two 
possible polarizations given by the equation 


r T 2 /2! L 
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Further, the index of refraction (n) is given by 


(3.7) 
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Equation (3.8) is known as the Appleton-Hartree formula. The above equations 
define a unique pair of values of polarization and index of refraction for 
each travelling wave mode. 

For the midlatitude launches from Wallops Island, Virginia, the earth's 
magnetic field is only 20° from a downward orientation (i.e., the -z 
direction). This implies that the propagating wave is nearly longitudinal 
with respect to the magnetic field. It follows that 



<< \l-X-jZ\ 


(3.9) 


Using this fact equations (3.7) and (3.8) can be simplified, yielding 


R • ± J 

2 _ . X 
n ‘ 1 -JZ±Y L 


(3.10) 

(3.11) 


This is known as the quasi-longitudinal approximation. A further approxi- 
mation can also be made in equation (3.11) when the local plasma frequency 
is much lower than the radio frequencies used in the propagation experiment, 
thus from equation (3.4) it follows that 
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Taking advantage of this fact, equation (3.11) becomes, by binomial expansion, 

„ . , h I 
k “ 1 “ 1 -jZ±Z L 

or, on rationalizing. 
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Equations (3.10) and (3.13) when used with the minus sign represent the 
parameters for the ionospheric mode called the ordinary mode. Similarly 
with the plus sign these equations represent the parameters for the second 
mode called the extraordinary mode. The ordinary mode has right-handed 
circular polarization while the extraordinary has left-handed. Equation 
(3.13) has the form of 


n R = "S ‘ 3 «/ 


n L * V L ~ 3 *1 


(for the ordinary mode) (3.14) 

(for the extraordinary mode) (3.15) 

When these are substituted in equation (3.1), the travelling waves become 
* . -Xpfcs j(o >t-v k 2 ) 

E R (* 3 t) = e R ° e R ° (3.16) 
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where k Q is the free space propagation constant. It is seen that x 
represents attenuation per unit length and u represents phase per unit 
length. 

Using the above analysis, Faraday rotation in the ionosphere can now be 
explained. To begin, assume that a linearly-polarized plane wave is 
propagating in the +s direction as in Figure 3.1. Also assume that it is 
incident on a semi-infinite ionospheric region of free electrons having a 
uniform magnetic field directed in the -z direction. Because the magnetic 
field is longitudinal to the direction of propagation, the travelling wave 
must be decomposed into the modes defined in equations (3.16) and (3.17) as 


& *V 
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it travels through the ionosphere. Any linearly-polarized wave may be 
directly decomposed into the vector sum of equal amplitude right- and left- 
handed polarized electric field vectors, which are the ordinary and extra- 
ordinary modes, respectively. Figure 3.2(a) shows this situation graphically 
where it is assumed the two polarized waves are initially in phase at z - 0. 
After the two modes have propagated a small distance. As, they are no longer 
in phase because their propagation constants are different. The effect of 
this is to cause the axis of polarization to rotate, which is the Faraday 
effect and is shown in Figure 3.2(b). The amount of rotation caused by 
this incremental change in position is 

k 

Afl = Y ( V y L )As (3 * 18) 


where a positive sign indicates rotation in a left-handed sense. 

For the purposes of the rocket radio propagation experiment the above 
analysis is not complete until the effect of the moving rocket is included. 
For a rocket moving in a +s direction, this can be done by expressing 
equation (3.18) in a differential form 


dt ' 2 


/ ^ dz 

V*! 5 it 


(3.19) 


Equation (3.19) can be rewritten in the following form 
1 


F -ii 'W 


-1 


where F is the Faraday rate in deg s , v is the rocket velocity in m s 
and 


(3.20) 
-1 


k R ” k o v R 

(3.21) 

k L ' k o v l 

(3.22) 


It follows from equations (3.4), (3.13), and (3.20) that Faraday rate is 
directly proportional to electron density. Therefore, experimentally 
measuring Faraday rate as a function of altitude can provide information 
on electron density also as a function of altitude. 

It should be pointed out here that although approximations have been 
made in the preceding analysis to allow relatively simple formulation of 
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Rpsultoni is still linearly 
polarised but iis plane of 
poor i/ation has been 
rotated through an angle Ail 


Figure 3.2 Illustration of the rotation of the plane of 
polarization of : he resultant wave as it 
propagates from a) s = 0 to (b) s = As along 
the magnetic field [Yeh and Liu , 1972]. 
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Faraday rotation, these approximations are not used in the actual experiment. 
3.3 Measurement Technique of the Original Propagation Experiment 

This section describes the Faraday rotation measurement used in the 
original radio propagation experiment. Also a general system and operational 
description is presented that will be useful later in understanding the 
details of the modified experiment. 

The original propagation experiment measures both Faraday rotation and 
differential absorption to obtain electron-density information [Mechtly et al.j 
1967] . To accomplish this the experiment uses a transmitted radio signal 
that easily allows both effects to be identified. This is done by indepen- 
dently generating two signals having the polarizations of the ionospheric 
ordinary and extraordinary modes and then combining them to form a composite 
transmitted signal. Figure 3.3 graphically describes the signals which are 
transmitted from a location near the launch site. This diagram represents 
what is seen at the rocket looking back towards the transmitting source. 

Figure 3.3(a) is the right-hand polarized (RHP) ordinary mode signal while 
Figure 3.3(b) is the left-hand polarized (LHP) extraordinary mode signal. 

It will be noted that the ordinary mode has a larger amplitude. The 
ordinary to extraordinary power ratio is adjusted to be 10 db; the reason 
for this will be discussed later in the section. Also, the extraordinary 
mode signal is adjusted to be nominally 500 Hz higher in frequency than the 
ordinary mode. 

When the electric fields from these two signals are added, the resultant 
shown in Figure 3.3(c) is obtained. This composite signal, seen at the 
rocket, is an electric field vector that spins in a right-handed sense at 
the mean of the ordinary and extraordinary frequencies. Because the 
amplitudes of the two circularly polarized signals are different, the 
composite signal is elliptically polarized. Further, since the two signals 
have different frequencies, the polarization ellipse spins in a left-handed 
sense at a rate of 250 Hz. As this composite signal propagates through the 
ionosphere, Faraday rotation will cause the phase of the polarization 
ellipse to increase with altitude according to the electron-density profile. 

In the reference plane of the moving rocket this change in phase appears as 
a change in frequency of the polarizati rotation rate. By measuring the 
frequency of polarization rotation and extracting the contribution due to 
Faraday effect the electron density can be determined. 
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Figure 3.3 The ordinary mode signal (a) and the extraordinary mode 
signal (b) combine to form the polarization ellipse fcl 

[Salah and Bowhill 3 1966]. ^ J 


The complete propagation experiment consists of two systems operating 
simultaneously, but on separate radio frequencies. The choice of operating 
frequencies is made to provide significant daytime differential absorption 
while still leaving sufficient signal for Faraday rotation measurements; 
the 2 to 5 MHz range has been found to adequately satisfy the above require- 
ments. Typically one system operates at 2.225 MHz, while the other operates 
at either 3.385 or 5.040 MHz. The criterion for choosing between the latter 
two frequencies is based upon which frequency combination is expected to 
give the best tradeoff between measurement sensitivity and highest: attainable 
measurement altitude. Measuremert sensitivity is inversely related to the 
square of the frequency, but lower frequencies have lower reflection heights 
beyond which the experiment ceases to operate. In general 3.385 MHz is used 
f night and 5.040 MHz in the day. 

The design of the original experiment was _ased on the best available 
technology of the early 1960’s. The approach taken was to keep the rocket 
payload simple at the expense of a more complicated ground transmitting 
system. The resulting design wes capable of supplying Faraday rotation 
resolution to ±5 deg s" 1 and absorption resolution to ±0.2 dB. When combined 
with Langmuir probe data, electron-density profile accuracy of better than 
±10% could be obtained [ Mechtly , 1974]. However, the above results could 
only be obtained during the daytime. 

The rocket payload equipment required to receive th** composite trans- 
mitted signal of Figure 3.3(c) includes a linearly polarized ferrite loop 
antenna feeding an AM-type receiver. That is, the receiver detector is of 
the envelope detection type and follows the carrier's envelope. For the 
original experiment this type of receiver works well with the composite 
signal because the 10 dB power ratio between ordinary and extraordinary 
modes (32% modulation) prevents the signal, as seen at the spinning 
rocket, from ever undergoing any deep modulation nulls. If a modulation 
null did occur, the envelope detector would cause intermodulation of the 
signal preventing proper operation of the experiment. 

The receiver output consists of DC and 500 Hz signals having an 
amplitude ratio corresponding to the 10 dB power ratio of the two modes. 

As seen in Figure 3.4, the output of the receiver is telemetered to the 
ground where it is recorded and also sent to the transmitter van location. 



CbL Von * * Woc*»t Payload 



Figure 3.4 Radio propagation experiment system block diagram 
[Knoebel et al. 3 1965], 
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At the van, the percent modulation of the receiver output is measured. 
Under conditions of differential absorption the extraordinary mode is 
absorbed more than the ordinary mode which causes the percent modulation to 
decrease below its nominal 32% value. When this is observed in the van, the 
extraordinary transmitter power is automatically increased to condensate for 
it. Power levels are recorded during the flight providing a measure of 
differential absorption. 

In addition to recording the transmitter power levels and the receiver 
outputs, the 500 Hz difference frequencies observed in the transmitter van 
and the rocket spin rate are also recorded. Faraday rate is measured from 
each system after the flight by measuring the 500 Hz difference frequency, 
fj, the rocket spin rate, f , and the receiver output frequency, f , and 

u S O 

using the relation 

F = h[f o~ f d~ 2f s ] (3 ’ 23) 

where F is the Faraday rate. 

Several techniques have been used through the years to make the 
necessary frequency measurements. Initially a mechanical system was used to 
directly obtain the Faraday rate [Gooch et al. , 1966] . Later a technique of 
counting zero crossings was used for frequency measurement [Slekys and 
Mechtly , 1970]. Recently, to give better discrimination against noise, a 
method of using spectral measurements from Fourier analysis has been used. 
This later approach is also used in the modified propagation experiment and 
is described in a later chapter. Further information on the original 
experiment can be found in Knoebel et al. [1965], Gooch et al. [1966], 
Knoebel and Skaperdas [1966], and Morrison [1966], 

3.4 Principle of Measurement Using Modified Propagation Experiment 

This section contains a general description of the modified radio 
propagation experiment; the details will be given in the following chapters. 
The purpose here is first to describe the signal that is transmitted from 
the ground installation to the rocket and to show how ionospheric effects 
modify this signal as received at the rocket. Then the modification of the 
signal by the receiver in the rocket payload will be discussed. 

3.4.1 Signal at the input of the receiver. The rocket, while it is 
ascending in the region of interest, has a velocity vector, v, that to a 
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good approximation lies along the line of sight from the ground transmitting 
installation; see Figure 3.5(a). The receiving antenna is sensitive to 
electric field components that are perpendicular to the rocket axis, which 
is usually nearly parallel to v. Figure 3.5(b) shows the field components, 
of the modified experiment, as they are observed from the rocket's location 
looking back towards the transmitter. The two linearly polarized CW signals 
emanating from the same dipole antenna are decomposed into right- and left- 
hand polarized electric field vectors; these are shown in the figure as 
and A £, respectively. The numerical subscripts differentiate between the two 
signals, which are at slightly different frequencies, f. and where 
(/j-fj) e 500 Hz as in the original experiment. Because the two CW sources are 
incoherent, their respective right- and left-hand mode vectors coincide along 
the polarization axis at different times. However, because both signals are 
transmitted from the same antenna and with the assumption that their frequen- 
cies are close enough, they both share the same polarization axis. As the 
rocket ascends through the ionosphere, this polarization axis will rotate due 
to Faraday rotation; in the northern hemisphere an observer positioned as shown 
in Figure 3.5(a) will see the polarization axis rotate in a left-hand sense as 
the rocket ascends. Finally, a phasor is included in Figure 3.5(b) that 
represents the rocket receiving antenna with effective length, l. This antenna 
rotates at the rocket spin rate, f , in a right-hand sense. It is only sensi- 
tive to that component from each of the four electric field vectors that is 
parallel to the antenna phasor at each instant of time. 

A mathematical description of the rocket receiver input voltage, V(t ) , 
can now be derived. To do this, the transmitted signal will be modified in 
stages to include the rocket's spin, forward velocity, and the ionospheric 
effects to be measured. This is begun by providing a mathematical description 
of the two independent CW signals as seen by a stationary observer positioned 
at some far-field altitude, z. 

As was stated previously each of the linearly polarized CW signals can be 
described as the sum of right- and left-hand polarized waves. In the 
immediate vicinity of the observer, these can be approximated as plane waves 
travelling with propagation constants, and k r , for right- and left-hand 
polarization, respectively. Mathematically, the electric fields of the two 
indepenuent CW signals can be described as: 
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Figure 3.5 Transmitted signal decomposed into right- and left-hand 

polarized components at frequencies and f2 as seen from 
the rocket. No differential absorption case ( A ^ = A L ) shown. 
f s is the roll rate of the rocket. 
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cosfoijt-fe^a+ep + y sinCojjt-fe^+ej)] 

♦ A R ^[x cosfiOjt-fe^a+ep - y sin^t-k^a+ej) ] 
and 

#2 = A r 2^ cos ^2 t ~ Ji R Z+e 2^ * $ sin^t-k^a+e^] 

+ A^ 2 [x cos^t-k^ 2 *^ " ^ sinfw^t-^s+e^) ] 

where 

Uj = 2n /j and 1^2 - 2irf 2 are the carrier frequencie: 

k R is the right-handed (ordinary) moue propagation constant 

kj- is the left-handed (extraordinary) mode propagation constant 

9j and ©2 are the carrier phases 

A R y and A r2 are ^e right-handed amplitudes 

A t , and A 7 are the left-handed amplitudes 


It was assumed above that since then k R ^ = k R 2 = k R and k^ 5 k ^ = kj. 
The error caused by this assumption is at least an order of magnitude smaller 
than the. desired accuracy. 

The total transmitted signal is 



+ 


The frequency spectrum of ^ is shown in Figure 3.6(a) where the amplitudes 
of the spectral components are shown. It can be seen that the ordinary and 

4 * 4 “ 

extraordinary modes for E ^ and E 2 are at the same frequency. 

The spinning rocket antenna can be described mathematically as 


A = l[x co s(u> t+<p) + y sin(w t+4>)] 
s s 


where 

l is the antenna effective length 
u) = 2irf is the rocket spin rate 

p is the rocket spin phase constant 




(b) Spectrum seen by 
rotating observer 
below the ionosphere 
(f = rotation rate 
in'ccw direction) 



Cc) 


Spectrum seen by 
rotating and linearly 
moving (away from 
source) observer 



(d) 


Same as (c) above, 
but with 0 and X 
experiencing different 
phase velocities, as 
in the ionosphere 



(e) Final spectra from 
original experiment 



Figure 3.6 The signals received at the rocket are frequency 
shifted due to: (b) rocket spin, (c) Doppler 

shift and (d) dispersive (or differential) 
Doppler shift. In the original experiment only 
the O 2 and components are transmitted. 
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It has been assumed above that the rocket antenna spins in a right-handed 
sense with respect to its direction of travel, which in this case is in a 
positive-s direction, at a spin rate f . 

The RF voltage at the input to the receiver is simply the projection of 
ffj, along ?, that is 

V(t ) = ^4 

Therefore 

V(t) = £{ cos (u> g £■*■<}> ) [^cosCwjt-fc^a+ej) + A^cosim^-k^+Q^) 

+ A R2 cos{u 2 t-k R z+Q 2 ) + A L2 cosiu 2 t-k L z+e 2 j] 

<■ sin((o g i+<(i)[4 i?1 sin(a) 1 ^-A: B 3+e i ) - /^sinfm^-^s+ep 

+ A R 2 sin{u 2 t-k R z+B 2 ) - ^sinfa^-k^s^)] } 

Using trigonemetric identities F(t) can be put into a form that shows its 
spectral components 

V(t) = Z{4 /?1 cos[(a) 1 -a) g )i-A: J? 2+0 1 -ji] + ^jCostfaj+w^iM^+e.^) 

+ A R2 cos[eu 2 -u> s )t-k R z+B 2 ~$) * 4£ 2 cos[(w 2 +u g )t-k£S+e 2 +i|>] } 

The frequency spectrum of vet) is shown in Figure 3.6(b). As can be seen, 

the effect of a right-handed (counter-clockwise) spinning observer is to 

pull the ordinary mode down in frequency and the extraordinary mode up in 

frequency by the spin rate, f . 

s 

To complete this mathematical description of the modified experiment, 
one must also include the effects of the rocket's motion. Typically, the 
rocket velocity, v, is between 1.0 and 1.5 km s - * at the altitudes of 
interest. There are two important effects that should be included. First, 
assuming that the rocket moves in a positive-a direction, it follows that z 
in the above equations, becomes a function of time, i.e., z t. Second, 
since the rocket is moving away from the transmitter, there will be a down- 
ward Doppler shift in the received frequencies. 

The relation between the transmitted frequencies and f 2 and the 
received Doppler-shifted frequencies, and f' 2> can be quantified. For an 
observer receding from a source, the equation relating / to f is 



' if*. •«'•**> t 
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r = /a - b 


where o is the speed of light. For typical rocket velocities and operating 
frequencies, the observed frequency is lowered by 10 to 20 Hz. This absolute 
shift in carrier frequency is of little consequence because the experiment 
is based on comparison of signals having nearly the same frequency which are 
all Doppler shifted. However, the differential Doppler shift between signal 
components may be of concern. The formula relating the difference frequency, 
f M - /“i - to its differential Doppler counterpart, = f[ ~ * s 

r M - V 1 - b 

F° r fu - 500 Hz and v between 1.0 and 1.5 km s *, f' would be between 
499.9983 aw 3 499.9975 Hz. If, as in the original experiment, it is assumed 
that fy = f’ the measured Faraday rate will be low by 1 to 3 mHz. This can 
cause an error of 10 to 100 cm ^ in electron concentration, which can be of 
significance in the nighttime ionosphere. The effect of Doppler shift on 
the signal spectra is shown in Figure 3.6(c); the small differential Doppler 
shift cannot be 'bserved on this scale. 

When the above effects are included in the analysis, V(t) becomes: 

7(t) = H^jCOsf + /l il cos[ (a)j+io g -fc z .y) t+e 1 +< t>] 

+ A /?2 cos[(u|-a) s -fc i? y)<:+e 2 -(|>] + 4 i2 cos[(a)'+u> g -fc L u)t+e 2 +<f>] } 

(3.24) 

The ionospheric phase velocity information, which will be used to retrieve 
electron-density information, is now contained in the frequency parameters 
of each of the sinusoidal arguments as k v and k v. This spectral modifi- 
cation is included in Figure 3.6(d); as is implied, the nighttime ionosphere 
provides only small changes in the observed spectra. 

When a similar mathematical analysis is done for the original experiment, 
its spectra are ideally as shown in Figure 3.6(e) where there is only one 
ordinary and one extraordinary spectral component present, having different 
amplitudes. It is interesting to point out that a practical limitation of 
the original experiment has been its inability to generate perfect circular 
polarization at each frequency. The result of this is that in practice, 

0 j and Z 2 components present in Figure 3.6(d) were also present, with reduced 
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amplitudes, in the spectra of the original experiment; at times this has 
caused problems in the experiment's operation. 

3.4.2 Signal at the output of the receiver. The preceding section 
presented a graphical and mathematical description of the signal in the 
modified experiment as it appears at the input of the receiver included in 
the rocket payloac. In this section the analysis is extended to the output 
of the receiver. It is this signal which is telemetered to the ground and 
Fourier analyzed to give the Faraday rotation rates from which electron 
concentration is deduced. 

The receiver's detector is of the envelope detection type. Envelope 
detectors are primarily designed for applications where the received signal 
is amplitude modulated (AM). For other modulation types, such as single- 
sideband (SSB) and double-sideband (DSB) , the envelope detector cannot be 
used without introducing non-linear intermodulation effects; for the case 
where the modulating signal is a complex waveform, such as voice, this 
intermodulation would have disastrous effects. Even so, if the received 
signal has a spectrum with discrete spectral components, as is the case for 
the modified radio propagation experiment, it is still possible to identify 
the spectral components at the output of the detector. This is because a 
discrete spectrum in turn generates another discrete spectrum, albeit more 
complicated, and providing adequate frequency resolving techniques are 
utilized, this more complicated spectrum can still be used for Faraday 
rotation analysis. 

To understand how the receiver input signal, derived in Section 3.4.1, 
is modified in the detection process, the envelope detector's operation must 
be explained. As the name implies, this type of detector is designed to 
follow the carrier envelope. The modulated carrier can be represented by 
the general expression 

V(t) = Kit) cos(wt) (3.25) 

The output of the envelope detector is proportional to the carrier's time 
varying amplitude, i.e., |#(£)|. 

It will be noted from the above, that the carrier envelope is such that 

' Kit), 
k -K(t), 


carrier envelope - 


Kit) > 0 
Kit) < 0 
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This implies that sign ambiguities can result, without additional information, 
when attempting to extract K{t ) from the carrier envelope, |X(t) | . This 
situation is illustrated in Figure 3.7 where the carrier envelopes for the 
particular examples of AM and DSB modulation are identified; the modulating 
source is a pure sine wave. As can be seen X(t) cannot be directly obtained 
from the DSB carrier envelope. However, because the envelope detector 
modifies K{t ) in a predictable fashion, it is still possible to obtain 
pertinent information about K{t) from the detector output, such as amplitude 
and frequency, since it is already known that K(t ) is a pure sinusoid. In a 
similar fashion, because the receiver input for the modified experiment is 
knowri, it is possible to identify its carrier envelope, which is the 
envelope detector output. In doing this, the spectral components at the 
detector output can be identified, which will allow the extraction of the 
needed Faraday rate information. 

To determine the carrier envelope of the receiver input signal equation 
(3.24) must be modified to have the form of expression (3.25). This can be 
done, in two steps, using trigonometric identities and the fact that 
a)' = First, equation (3.24) is modified as follows 

V(t) = i{[A i?1 cos0 fll+ 4 z . 1 cos0 L1+ 4 i?2 cose i?2+ A i2 cose i2 ]cos(a)'*) 

(3.26) 

+ ^l sine flr 4 Ll sin0 Ll + ^2 sin0 7?2 + ^2 Sine L2 ]sin(a) l t)} 

where 

e ifi * (w^-V* 

e Ll ' 

e s2 * t“i + VV )t ‘Y* 

e i2 ■ (“m-vV'-V* 

Second, two new functions, k(t) and a(t), are defined as follows 

*(i)cosa(t) = ^iCos0 i?1 ^ L1 cos0 2 . 1 ^ i?2 cos0 i?2+ 4. z . 2 cos0 L2 

*(*)sina(t) = /l /fl sin0 /?1 -yl j: . 1 sin0 L1+ 4 /?2 sin0 i?2+ 4 £2 cos0 i . 2 


(3.27) 

(3.28) 
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If equation (3.27) is substituted into equation (3.26), the result is 


V{t ) = 7fe(£) cos [u> '£-<*(£)] 


(3.29) 


As can be seen, except for a multiplicative constant that is not critical to 
the analysis, the carrier envelope can now be identified as |#(t) | . The 
phase modulation term, a (t) , is of no consequence to the envelope detector 
which is not sensitive to phase or frequency modulation effects about the 
carrier frequency. 

From this analysis a specific equation for the detector output can be 
derived. Combining equations (3.27) and (3.28) it follows that 


X 2 (t) « (fr(t)cosa(t)] 2 ♦ [Ji(t)sinaft)] 2 


^fll * A Ll * A R2 * A L2 1 * 2J ffl' 1 Il cos t e ffl* e £p 


' U Rl A S2 cas( ^2-W * U 

* M lA2 cos(e i,l* e i!2> * 24 

* 2 ‘ l S2 /1 £2 cos(8 i;2- 8 I2 ; 

By further algebraic manipulation we have 


R1 A L2 

L1 A L2 


C0S(6 L2 V 

C0S ^L1 +Q L2 ) 


(3.30) 


K 2 m = lA R1 *t L1 >t B2 *A L2 } 2 * 

+ 24 /?l /1 i?2^ COs(e i?2" 0 i?i :i “ 1 ^ + 2A Rl A L2^ COS( ' Q L2' & Rl ) ~ 1 ^ 

+ 2A Ll A R:^ COs( - d Ll* B R2^~ 1 ^ + 2i4 Ll' 4 L2^ COs( ' 0 Zil +e L2- ) " 1 'l 

+ 2A R 2 A L2^ COS ^ Q R2~ Q L2^ ( 3 . 31 ) 

2 

The next step is to determine the square root of K (t) . An appropriate 

2 

method to do this is by using a binomial expansion of K (£). The particular 
2 

form of X (t) , shown in equation (3.31), was chosen because it satisfies the 

necessary convergence conditions of the square-root binomial expansion. It 

is only necessary to carry out the expansion to the second degree 
h 2 

([1+X]^ = 1 + X/2 - X /8 + ...) to demonstrate and identify the relevant 
features of the detector output spectrum. However, because of the length of 


32 


even this much expansion, the final analytic expression for the detector 
output will not be given here; instead, a partial listing of the output’s 
spectral components, resulting from the second degree expansion, is given 
in Table 3.1. For each frequency listed, its corresponding phase and 
approximate amplitude is given. The amplitude information is only approximate 
because each degree of the expansion has terms that contribute to these 
amplitudes; therefore, only the dominant amplitude term is shown for each 
spectral component in the tabulation. 

Table 3.1 shows that the spectral components with the largest amplitudes 
are and The components ?t multiples of [2/ g + (l/2ir ){k R -k^)v] away in 
frequency from these have decreasing amplitudes. Figure 3.8 graphically 
•presents an actual receiver output spectrum taken from early in the flight 
of Nike Apache 14.533. As can be seen, the predicted spectrum is quite 
evident . 
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4. MODIFIED PROPAGATION EXPERIMENT 

This chapter is divided into three sections, corresponding to the 
three major parts of the experimental system. First/ the ground-based 
transmitter system will be described. Second, the experiment payload 
package will be presented. Finally, the telemetry and data recording 
systems will be briefly explained. 

4.1 Ground-Based. Transmitter 

The transmission source for the dual CW signal described ir. the 
previous chapter will be detailed in this section. As in the original 
experiment, two systems are simultaneously operated at different 
frequencies; system 1 operates at 2.225 MHz while system 2 operates at 
either 3.385 or 5.040 MHz. The two systems are not identical, as they were 
developed for the original experiment at different times. Therefore, there 
is some difference in system configuration between systems 1 and 2 for the 
modified experiment; pertinent differences will be identified in the 
following discussion. 

Block diagrams for systems 1 and 2 are shown in Figures 4.1 and 4.2, 
respectively. When compared with the block diagram of the original 
experiment (Figure 3.4), the reduction in system complexity is apparent. 

The fundamental difference between systems 1 and 2 is the method used to 
generate the independent carrier signals, at frequencies and Once 

these signals are generated and combined, both systems operate identically 
by providing power amplification and transmission using a half-wave dipole. 

The technique used in system 1 to generate the carrier signals is to 
use two separate crystal oscillator stages. Both crystals are housed in 
the same linear-proportional crystal oven. Placement of the crystals in 
the same oven helps reduce relative drift between the two frequencies. 

The use of a linear-proportional oven is preferred to an on-off type, as 
the latter causes sudden jumps in “he frequency difference, fy = 
whenever the oven is switched on or off. The output of each oscillator is 
used to drive a buffer-amplifier stage before being combined. Additionally, 
the oscillator outputs are sampled and then mixed together to form the 
500 Hz reference difference frequency, fy. This reference frequency is 
both monitored at the transmitter and sent, on data lines, to the main 
base telemetry station where it is recorded on analog and digital magnetic 
tapes for later analysis. 




Figure 4.1 System 1 configuration. 
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System 2 follows a slightly different approach to generate and 
As shown in Figure 4.2 there are again two independent crystal oscillator 
stages whose crystals are located together in a crystal oven; their 
frequencies are 1000.250 kHz and 999.750 kHz. The outputs are directly 
added together, using a power combiner. As in system 1, the difference 
frequency is also generated for monitoring and recording purposes. The 
pov. ex -combined signal is then balance-mixed with another local oscillator. 
The local oscillator frequency, / , is chosen so that the sum frequency, 

/ +1 MHz, is either at 3.385 or 5.040 MHz, whichever is desired. The output 
of the balance-mixer goes to a buffer-amplifier stage, which in turn is 
used to drive the linear power anqslifier. 

Up to this point, the value of the reference difference frequency, f 
has been quoted as 500 Hz. However, for the modified propagation experiment 
any frequency between 400 and 530 Hz would be acceptable, as it is 
not critical to the experiment that is exactly 500 Hz. It has been found 
to be critical though, that the difference frequencies from systems 1 and 2 
should not be the same. The reason for this is because there is cross-talk 
interference in either the data lines or equipment at the main base 
teiemetry station prior to recording the reference frequencies for each 
system. If the frequencies are nearly the same, they can introduce error in 
the post-flight measurement of these reference frequencies for data analysis 
purposes. It has been established that there should be a minimum 35 Hz 
separation between the two reference frequencies to provide adequate 
protection against cross-talk on any of the recorded data. 

The power combiner, used in both systems is a Merrimac Industries power 
splitter/ combiner. It is a broad-band device specified to cover 2 to 32 MHz, 
although tests demonstrated that it worked satisfactorily at the 1 MHz 
frequency used in system 2. These units were also specified to have a 1 W 
maximum input power; they were found to also operate satisfactorily at a 
2 W power level in system 1. The signal level was adjusted from each 
oscillator sc that equal signal levels were observed, on an oscilloscope, 
at the output of the combiner. 

The original propagation experiment used two linear amplifiers in each 
system to provide the necessary transmitted power level; the modified 
experiment uses only one amplifier in each system. These amplifiers are 
specified to deliver 1 kW for an input of luO mW; in practice, however, 
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the maximum output is about 500 W. Each amplifier has an adjustable 
automatic level control, which is set at its minimum sensitivity position 
to prevent any amplifier- induced distortion. The input level to each 
amplifier is controlled by a Kay step attenuator and it is adjusted for 
maximum output power without overdriving the amplifier. 

The output of the amplifier from each system is used to feed a half-wave 
dipole. The dipole for each system is placed between telephone poles 
approximately 90 feet high and is fed with RG-8 coaxial cable. A balun 
(matching transformer) is used at the center of the dipole to help prevent 
radiation from the coaxial cable which would distort the far-field antenna 
pattern. The antenna input impedance as seen through the coaxial feedline 
is typically 50 to 60 ft at an angle of ±20°. 

4.2 Payload Instrumentation 

The rocket payload used in the modified radio propagation experiment 
is identical to that used in the original experiment. The payload section, 
used in the propagation experiment, will now be described. 

The primary payload section, related to the propagation experiment, 
consists of two receivers used to receive the ground-based transmitted 
signals from systems 1 and 2; the receiver in current use was developed at 
the University of Illinois Aeronomy Laboratory by G. W. Henry, Jr. A 
photograph of three completed receivers is shown in Figure 4.3; note that 
one has its covers removed to show the electronic circuits. These 
receivers were designed to provide good selectivity with a linear detector 
and AGC response over a wide dynamic range in a rocket environment. 

Receiver specifications are shown in Table 4.1 and its detector characteris- 
tics are shown in Figure 4.4. 

The antennas used by the payload receivers are magnetic dipoles, and 
their design is typical of that used in transitorized portable broadcast- 
band radios. Each antenna consists of two ferrite rods used in the 
core of a Litz wire coil which is in parallel resonance with variable 
capacitors. A shielded one-turn symmetrical loop is used to couple the 
antenna to the receiver input. A schematic of the antenna circuit is shown 
in Figure 4.5. The payload section, housing the two antennas is shown in 
Figure 4.6. It will be noted from this photograph that each antenna is 
placed to maintain symmetry between the two ferrite rods; this is done to 
keep the antenna pattern as symmetrical as possible. 
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Table 4.1 Receiver performance data 


Input : Frequency 

Level 
Impedance 


Selectivity: -6 dB bandwidth 

-60 dB bandwidth 
Shape factor (60/6) 

Spurious responses (for 2.225 MHz f c ) : 
(2.0 to 5.0 MHz range) 


AGC Characteristics: 

Detector Output Flatness: 
AGC Voltage Range: 

AGC Output Impedance 

Detector Characteristics: 

Nominal dc Output 
Linearity 

Temperature Stability 
(-10 to +50°C) 

Output Impedance 

Power Requirements: 


Environment: Temperature 

Shock and Vibration 


2.0 to 5.0 MHz 
-110 to +10 dBm 
To match ferrite antenna 
(600 0 nominal) 


2.811 

kHz 



4.490 

i kHz 



1.597 

: 1 



2122 

kHz 

-103 

dB 

2156 

kHz 

-103 

dB 

2451 

kHz 

-102 

dB 

2470 

kHz 

-100 

dB 

3135 

kHz 

- 50 

dB (Image) 


±0.2 V, -110 to +10 dBm input 
+5.0 V (minimum signal) 

0.8 V (maximum signal) 

2.2 kO 


2.0 ±0.2 V 
±2 dB, 0.15 to 4.0 V 
Gain: ±3 dB 

Offset: ±0.05 V 

Less than 10 ft 

+20 to +35 V @ 55 mA 
-20 to -35 V 0 15 mA 

-10 to < S0°C 
0 to 15 g, sinusoidal 
30g , random 
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In addition to these payload receivers, the propagation experiment 
makes use of a magnetometer to measure the rocket spin rate. The magnetom- 
eter is an instrument sensitive to the earth's magnetic field. Using this 
field as a reference, the magnetometer has a sine wave output with frequency 
equal to the rocket spin rate. 

Both receiver outputs and the magnetometer signal are telemetered to the 
ground where they are recorded for later analysis. The telemetry and data 
recording system will be described in the following section. 

4.3 Telemetry and Data Recording Systems 

For the purposes of data analysis, it is required that the two rocket 
receive?.* outputs, the magnetometer signal, and the two reference frequencies 
from the transmitters be recorded throughout the duration of the rocket 
flight. Therefore a telemetry and data recording system must support the 
experiment; details of both these systems are presented in this section. 

The telemetry data from the rocket is transmitted using a frequency- 
modulated (FM) carrier system where the modulating signal is a composite 
of all the data channels. To form this composite modulating signal, each 
data channel first is used to FM a voltage-controlled oscillator (VCO) . 

Each VCO operates at a different center frequency with enough separation 
so that the FM spectra from any one VCO does not interfere with the others. 
The VCO center frequencies are chosen from twenty-one standardized frequency 
channels (IRIG subcarrier channels: 1 to 21). Higher numbered channels 

require proportionally more bandwidth, but can carry higher frequency in- 
formation. Therefore, the particular channel number chosen for a data 
channel is determined by its highest frequency component to be transmitted. 
The two rocket receiver outputs use IRIG channels 16 and 17 which can pass 
frequencies up to 600 and 790 Hz, respectively. The 5 to 10 Hz rocket spin 
(magnetometer) signal used IRIG channel 6, that passes frequencies up to 25 
Hz. Once all the FM subcarriers are generated, they are simply linearly 
combined to form a composite modulating signal. The RF carrier frequency 
used at NASA-Wallops Island is in the vicinity of 2.25 GHz (S-band) . The 
transmitted power output is typically 2 to 3 W. 

Once received, the rocket telemetry signal is FM detected to obtain the 
composite modulating signal. An analog magnetic tape recording is made of 
this composite signal. Also recorded, on different tracks of the same tape, 
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are various multiplex signals containing time code information and the 
reference frequency signals from the transmitters. Additionally, a 
precision 100 kHz reference oscillator signal is recorded on the analog 
tape. This reference signal is used to compensate for the variations in 
tape speed during playback of the tape. Recording typically begins two 
minutes prior to launch and continues until splashdown. 

The analog tape is not generally used in the processing of propagation 
experiment data. Instead a digital tape recording is made in real-time 
that records digital samples from the five previously mentioned data 
channels 

a. system 1 receiver output (IRIG channel 17) 

b. system 1 reference signal 

c. system 2 receiver output (IRIG channel 16) 

d. system 2 reference signal 

e. magnetometer signal (IRIG channel 6) 

Since the signals from the two rocket receivers and the magnetometer are 
contained in the detected composite telemetry signal , they must be extracted 
using frequency discriminators tuned to the proper IRIG channels. The 
discriminator outputs have low-pass 600 Hz filters for the rocket receivers 
and a 25 Hz filter for the magnetometer. 

All five of the above data channels, a through e, are sampled 
sequentially. A complete set or sequence of samples is defined as one 
data frame. Five thousand data frames are recorded each second, which is 
equivalent to 5000 samples per data channel per second. The analog-to- 
digital converter has a resolution of 12 bits per sample. 

The tape format is such that the digital information on the tape is 
divided into groups (called records). There are three different types 
of records: header records, calibration records, and data records. The 

tape begins with 4 header records which provide rocket identification 
information; each header record contains not more than 45 words, i.e., 45 
digital values. This is followed by five telemetry discriminator cali- 
bration records, each containing 1005 words. Data records begin after the 
header and calibration records. Each data record contains 2008 words, of 
which 2000 are data samples and they are located in words 6 to 2005. Words 
2006 through 2008 contain time code information, in Universal Time, to the 
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nearest tenth of a millisecond that corresponds to the time of the first 
sample of the following data record; see Table 4.2. 

Occasionally a real-time digital tape is not possible. When this 
occurs, a post-flight digital tape is made from the analog tape. To 
eliminate the effects of tape speed variation from the analog tape playback, 
the sampling period is derived from the recorded 100 kHz reference signal. 
Also on occasion, the above described format is not achieved, therefore it 
is recommended that proper tape organization be confirmed prior to the 
tape's use. This can be easily done using a digical tape reading program 
developed by K. L. Fries, of the Aeronomy Laboratory Rocket Program. This 
program is documented in Appendix I. 



Table 4.2 Coding of time words [Fillinger al. , 1976] 


Word 2006 bits 

13 - 16 

9-12 
5-8 
1-4 

Word 2007 bits 

13 - 15 

9-12 
5-7 
1 - 4 

Word 2008 bits 

15 - 16 

11 - 14 

7-10 
5-6 


hundreds of milliseconds 
tens of milliseconds 
ones of milliseconds 
tenths of milliseconds 

tens of minutes 
ones of minutes 
tens of seconds 
ones of seconds 

hundreds of days 
tens of days 
ones of days 
tens of hours 


1 


4 


ones of hours 
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5. DATA REDUCTION 


The data collected from the radio propagation experiment is used to 
calibrate the Langmuir probe mentioned in Chapter 1. To accomplish this, 
several data reduction steps must be performed in sequence. These steps 
include the extraction of Faraday rates, the conversion of Faraday rate 
information to electron concentration, and the calibration of probe current. 
This chapter discusses each of these steps. 

5.1 Foiwier Spectral Analysis 

In Section 3.4 it was shown that the desired Faraday rotation information 
is embedded in the frequency arguments of certain sinusoidal terms from the 
rocket receiver output. Further, it was '•hown that the receiver output 
contained a large number of discrete spectral components with an approximate 
frequency spacing of 2f . Therefore, to extract the Faraday information, 

O 

a technique is required that can distinguish and measure the frequencies of 
particular spectral components. 

From Table 3.1, the receiver output frequencies of most interest arc / 
and f t where 


4 -- V - -4 - f. w 


4 * V * -4 * h w 


(5.1) 


(5.2) 


This is because the amplitudes of these spectral components are the largest 
of all the frequencies containing ionospheric information. Equations (5.1) 
and (5.2) each contain Faraday rate information as well as the Doppler-shifted 
reference frequency and the rocket spin rate. It will be noted that by 
subtracting from / the reference frequency is eliminated leaving only the 
rocket spin rate and the Faraday rate. The spin rate is sepai'atcly obtained 
from the magnetometer channel. Using the above frequencies, an expression 
can be given to determine the desired Faraday rate, F. 


p - i u\-r,-vj 


(5.3) 


FilLhun'v ct aL [J976] have described a procedure for deriving the 
frequency information from spectra of the signals. To understand how this is 
done consider an analog signal, x(^), having a general spectrum defined by 
its continuous Fourier transform (CFT) 



so 


fOO 4 

*(/) = X(t)e" j2ir ^ dt (5.4) 

J .00 

This definition, however, is not practical for the present situation because 
it provides no information on how the frequencies of equation (5.3) vary with 
time. A more useful definition is available for situations where the signal 
is quasi-stationary, i.e., for a signal composed of frequency terms that vary 
slowly with time, which is the case for the receiver signal. For the above 
condition, a short-term spectrum, X[f,t) can be defined as 

«ae> # 

X(f,t) = w(t-x) x( T )e' j2lT ^ T dx (5.5) 

J _ao 

where w(t) is a weighting or window function that has non-zero value for only 
a finite duration near t = 0 [Allen, 1977]; more will be said about this 
window function later. This new spectral definition allows for the possibility 
of observing spectral frequencies as they vary with time. 

Equation (5.5) has two equivalent but distinct interpretations [Allen 
and Rabiner, 1977]. One view is shown graphically in Figure 5.1(a) where 
(5.5) has been decomposed into two cascaded sections. The first section is a 
complex mixing operation, which has the effect of shifting the spectrum of 
x(i) down by an amount f. The second section is a linear filter witn impulse 
response, w(t). In general, W[f) , the CFT of w(£), has a low pass character- 
istic with some cutoff frequency, f . The output from the combined system is 
a complex time varying signal that is sensitive only to the spectral energy 
from the input that is present in the frequency interval (/-/ . /♦/* ). t’y 
altering the frequency parameter, f, the entire short-term spectrum can be 
decermined. 

The second interpretation of equation (5.5) is that it represents the 
CFT of the new signal 

y(x) = tf(t-x)x(x) (5.6) 

where t is now a fixed parameter. The transform of y[ x) is related to the 
transforms of x and u by the complex convolution of the two. 

y(f) = [e~j 2vft W[-f)]*X(f) (S.7) 

This point of view demonstrates that the short-term transform X(f,t) = Y{f) 
is greatly influenced by the choice of window function. This interpretation 
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LOW PASS FILTER 
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(a) 



Figure 5.1 Two interpretations of the short-term spectrum defined in 
equation (5.5). (a) System interpretation using a mixer 

and linear filter, where w(£) is assumed to be low pass. 

(b) Short-term spectrum is interpreted as the normal 
Fourier transform of x(t) multiplied by the sliding window, 
w(£-t) . 
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is shown graphically in Figure 5.1(b) where an arbitrary sliding window 
function is shown superimposed on x(t) ; the product of these two is then 
Fourier transformed. 

When the signal to be analyzed has been sampled to form discrete digital 
values, an alternative to the CFT must be used, this is known as the discrete 
Fourier transform (DFT) . The DFT is defined as 

00 

X(e 3W ) = l x(nT)e~ jm (5.8) 

n*-°° 

where x(nT ) is the sampled signal with sampling period, T [Oppenheim and 
Schafer, 1975]. The w term in equation (5.8) is the standard convention for 
the radian frequency transform variable which has dimension of radians. 

This should not be confused with the sinusoidal angular frequency term 
having dimension of radians per second. The relation between radian frequency 
and the corresponding sinusoidal frequency, f, is 

0 ) = 2 vf? (5.9) 

Also note that because the complex exponential is periodic, so also is the 
DFT defined by equation (5.8). In a similar fashion, a discrete equivalent 
to the short term spectrum can also be defined. 

. “ • 

X (e J “) = l v[(n-m)T]x(m)e~ 0um (5.10) 

n m=-°° 

Here again it can be seen that the window, uCnT ) , determines what portion of 
x(nT) is to be spectrally analyzed. Once a window function is chosen, 
equation (5.10) can be used to obtain the necessary spectral information from 
the receiver and magnetometer channels. 

Direct application of equation (5.10) to measure a signal's entire 
spectrum would require that it be repeatedly used, as each computation of 
(5.10) returns information for only one spectral frequency. To avoid this 
complication the actual procedure used for the propagation experiment data 
analysis uses a fast Fourier transform (FFT) to perform the DFT in equation 
(5.10) •[ Bergland , 1969]. The FFT is an efficient algorithm for computing the 
DFT that takes advantage of certain properties possessed by complex exponen- 
tials. By using the FFT there is a substantial reduction in processing time 
below what would be needed for direct DFT computation. The FFT uses an even 
number, N, of signal samples selected by the window function to calculate 
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N uniformly spaced components of the spectrum defined in equation (5.10); 
these components are at integer multiples [-#/2, -1, 0, 1, .... (N/ 2 -1)} 

of u = 2i/N. This corresponds to sampling the short-term spectrum at 
multiples of / = for example, sampling at a 5 kHz rate and using 5000 

samples would mean sampling the spectrum at 1 Hz intervals. 

Proper operation of the short-term spectrum analysis, using the FFT 
algorithm, requires that proper consideration be given to the data sampling 
rate, digital quantization effects, and the choice of window function. To 
satisfy the Nyquist sampling rate criterion, discrete samples must be taken 
at a rate that is at least twice as high as the highest frequency present in 
the spectrum of the signal being sampled. This implies that since the 
receiver channels have 600 Hz low pass filters, the sampling rate should be 
at least 1200 Hz. In practice since the low pass filters do not completely 
attenuate all frequencies above 600 Hz, a guard band must be provided, 
therefore, the experiment uses a sampling rate of 5000 Hz. If a signal is 
sampled below its Nyquist rate, a distorting effect (aliasing) results and 
causes high spectral frequencies to appear as low spectral frequencies. 

Some consideration should be given to the accuracy and resolution of the 
digital samples. The act of digitizing the data inherently introduces noise 
into the data and is known as quantization noise. The question of accuracy 
is less critical as absolute amplitude information is not necessary to the 
Fourier analysis. However, the analog-to-digital conversion should not 
introduce any non-linear compressional effects that might alter the signal 
spectrum. The present digitizing system used at the Wallops Flight 
Center (NASA) has a 12-bit resolution. 

The choice of window function is very critical to the success of the 
short term spectrum analysis. This can be noted by the predominant role 
that the window function plays in the short-term spectrum d ‘inition, both 
continuous and discrete. Although the particular window chosen has several 
effects [Allen, 1977] the predominant consideration for spectral analysis 
is the window's ability to reduce spectral leakage. Spectral analysis of 
finite extent data with duration NT requires that a projection of the data 
be made on a basis set spanning the observation interval. For Fourier 
analysis this basis set consists of trigonometric functions that have periods 
which are submultiples of NT. If the signal being analyzed contains only 
frequencies from the basis set, they will project onto individual basis 
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vectors. However, if the signal contains frequencies other than those from 
the basis set, then the signal is not periodic in the interval NT and will 
thus have discontinuities at the ends of the interval. These discontinuities 
are cause for the signal to project over the entire basis set; this is 
referred to as spectral leakage. The purpose, therefore, of the window 
function is to reduce these discontinuities and thus spectral leakage so 
that a signal will have significant contribution on only those basis vectors 
that are near the frequency of the signal. This is done by choosing a 
window that brings the signal to zero as smoothly as possible or, conversely, 
that has a transform with a good low-pass characteristic and low side-lobes. 

A large number of possible windows have been described in the literature, 
each satisfying the above criterion to varying degrees. An excellent survey 
and study of window functions for harmonic analysis using the DFT is given 
by Harris [1978]. In general it can be stated that the shape of the window 
function will determine the peak side-lobe level while the window duration 
will determine the transition length between passband and stopband. 

The actual window function used for the spectral analysis of the propaga- 
tion experiment data is the Gaussian window [Fillinger et at ., 1976]. This 
window shows the desired low-pass frequency response of a concentrated main 
lobe at dc with low level side-lobes at higher frequencies. Theoretically, 
the Fourier transform of a Gaussian function is another Gaussian and, there- 
fore, has no side-lobes in the frequency domain. However, because the FFT 
uses finite extent data, the Gaussian window function must be truncated at 
some point. This truncation causes the Fourier transform of the Gaussian to 
be slightly modified introducing the side- lobes. This situation is shown in 
Figure 5.2 where a 50 sample Gaussian window is shown with its DFT. The 
actual Gaussian window used for the data analysis has a duration of 5000 
samples and combined with the 5 kHz sampling rate, allows the FFT to sample 
the spectrum at 1 Hz intervals. 

The frequencies required for equation (5.3) must be known more accurately 
than are provided by the 1 Hz sample interval of the FFT. E. K. Walton has 
developed an algorithm to accurately interpolate between spectrum samples to 
determine the frequency of the spectral peak [Fillinger et al., 1976], His 
algorithm uses the fact that the Gaussian window modifies the spectrum of any 
signal having discrete frequencies or having a spectrum that is already 
Gaussian, in a known and predictable fashion. By using three spectral samples 
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Figure 5.2 (a) Fifty sample Gaussian window, (b) Log- 

magnitude of its corresponding discrete 
Fourier transform [Harris, 1978] . 
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about the peak, this algorithm is able to return the spectral peak frequency 
to within an experimentally determined accuracy of ±2.5 mHz. 

A computer program incorporating the FFT with Gaussian window and 
Walton's algorithm is used to extract the frequency information needed for 
equation (5.3). This program is documented in Appendix II. The program is 
used to separately measure both frequencies (/ and f ) from the receiver 
channels as well as the rocket roll rate (/ ) from the magnetometer signal. 

As examples of the application of the analysis Figures 5.3 and 5.4 show the 
measured Faraday rate versus seconds after launch for Nike Apaches 
14.533 and 14.534 determined from both system 1 (2.225 MHz) and system 2 
(3.385 MHz). 

5.2 Calculation of Electron Density 

Once the Faraday rate is obtained using the analysis of the previous 
section, the corresponding electron concentration value may be calculated. 

To accomplish this, an iterative computer analysis procedure is used 
[Mechtly et at . , 1970], This procedure is based on the complete Sen-Wyller 
theory [Sen and Wyller , 1960], which is an extension of the Appleton-Hartree 
analysis from Section 3.1. The principal change is the way that the velocity 
distribution of the electrons is treated. In practice the difference between 
the two theories is unimportant in the ionosphere at altitudes above about 
80 km. The Sen-Wyller theory is used here because the computer program was 
developed originally for use in the daytime ionosphere. 

Figure 5.5 shows a block diagram of the computer analysis procedure. 
Initially an assumption of electron concentration is made. Then using the 
Sen-Wyller equations and Burke-Hara polynomials [Mechtly et al. 3 1967] the 
refractive index is calculated for each ionospheric mode at the frequency of 
the experiment. With these indices of refraction and the rocket velocity 
the corresponding Faraday rate is calculated. This calculated rate is compared 
with the actual measured Faraday rate to determine a better assumed electron 
concentration value by using the equation from the final block of Figure 5.5. 
With the new value of electron concentration, the above procedure is repeated. 
This iteration is repeated up to six times to determine the actual electron 
concentration value corresponding to the measured Faraday rate. This analysis 
makes the assumptions that the direction of radio propagation and the rocket's 
direction of travel are along the line of sight between the launch site and 
the rocket. The first of these assumptions is in danger of failing when the 







Figure 5.4 Faraday rate for Nike Apache 14.534. A telemetry 

malfunction occurred between 60 and 68 seconds after 
launch. 




Figure 5.5 Procedure for analyzing Faraday rotation rates 
[Fillinger et at. , 1976] . 
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rocket is near the radio wave reflection height as there may be significant 
bending of the propagation path at this altitude, lhe second assumption 
fails near rocket apogee. 

When no differential absorption is measured, as at night, a collision 
frequency model must also be supplied. This model makes the assumption that 
collision frequency is proportional to the neutral atmosphere pressure, p. 
That is 

v = K • p 

where K is the proportionality constant with units of s” 1 Pa" 1 and is 

S -1 -1 

assumed to have a value of 6.3 x 10 s Pa [Smith et at 1978]. The 
atmospheric pressures are obtained from the COSPAR International Reference 
Atmosphere Model [CIBA, 1972]. The program that calculates the electron 
concentration is documented in Appendix III. 

The electron concentration values versus altitudes that were calculated 
for the flights of Nike Apaches 14.533 and 14.534 are shown in Figures 5.6 
and 5.7, respectively. Shown in each figure are the independent results 
obtained from the 2 MHz system and the 3 MHz system. Also shown in each 
figure, for clarity, is a partial plot of the calibrated pr' ! ; current. This 
calibration procedure is presented in the following section. 

5 . 3 Calibration of Probe Current 

The current from the Langmuir probe follows the detailed variation of 
the electron-density profile [ Meohtly et al ., 1967]. However, because no 
accurate equations exist that relate probe current I to electron concen- 
tration, N , the probe data must be calibrated by comparison with the results 
from the propagation experiment. To accomplish this, the ratio N g /I AVG is 
calculated for each electron concentration value. A listing of this computer 
program is presented in Appendix IV. 7he required inputs include up to sixty 
seconds of electron concentration values and the corresponding probe current 
data. The probe current must start ore second before the time of the first 
electron concentration value and must have all probe sweeps that occur above 
130 km removed. 

A semi-log plot is next made of ^ versus altitude. This is shown 

in Figures 5.8 and 5.9 for Nike Apaches 14.533 and 14.534, respectively. In 
general need not be a constant, however, any change is expected to be 

slowly varying with altitude. Therefore, an average ® e /I A y G value is 
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Figure 5.6 Electron-density values from the propagation experiments 
(individual points) compared with the calibrated probe 
current for i*ike Apache 14.533. 
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Figure 5.8 The probe calibration for Nike Apache 14.533. 
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Figure 5.9 The probe calibration for Nike Apache 14.534. 
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Figure 5.11 Electron concentration profile obtained by combining data 
from the probe and propagation experiments on Nike Apache 
14.534. 
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6. EVALUATION OF SYSTEM PERFORMANCE 

6 . 1 Introduction 

The purpose of this study is to identify sources of experimental error 
that presently limit the accuracy of the propagation experiment. The main 
discussion refers specifically to the modified version of the propagation 
experiment. However, much of the discussion is also relevant to the 
original experiment. 

The experiment, in both versions, is fundamentally concerned with the 
measurement of two receiver frequencies and the rocket roll rate included in 
equation (5.3), thus any sources of noise or anomalous behavior that affect 
these measurements are of importance. The results of this study are divided 
into those errors related to the hardware and those connected with the 
software frequency measurement procedures. 

6.2 Hardware 

6.2.1 No-u 2 sources . Noise is an inherent component of any electrical 
signal ;.nd can be loosely defined as any unwanted random or periodic effect 
that distorts the desired response. Noise of some form and magnitude is 
generated in all electrical components, including those that make up the 
propagation experiment. Further, radio frequency noise from a variety of 
naturally occurring sources is present in the reception of all free space 
radio propagation. It is, therefore, of interest here to identify any noise 
sources that are of sufficient amplitude to affect the accurac 1 of the 
experiment. Towards this end a series of system checks have been carried out 
at Wallops Island. These checks included an evaluation of the dependence of 
experimental accuracy on signal -to-noise ratio (SNR) and checks of noise 
effects which frequency modulate the reference signal. The checks have been 
made at various locations along the experimental system. 

6.2.2 Si jnal-to^noise ratio . Intuitively it may be expected that 
experimental accuracy must be dependent on SNR. In particular the relative 
level of signal compared with noise as observed at the telemetry station is 
of interest since this is where analog data is recorded for later processing. 
Possible noise sources that can contribute to the total noise level include 
those in the ground transmitting equipment, atmospheric noise, receiver noise, 
noise in the telemetry equipment, as well as noise in the data recording 
process. To test the dependence of accuracy on SNR, the propagation experi- 
ment system was set up to simulate pre-launch conditions. The radio trans- 
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mitter was adjusted for the typical power of 250 W, the rocket was mounted 
on the launcher in proper launching position and the telemetry signal was 
transmitted to the telemetry station. With this configuration digital and 
analog data recordings were made from the received telemetry signal for periods 
of 30 seconds using the procedure of Section 4.3. At the end of each period 
of 30 seconds the transmitted power was reduced by 10 dB and another 30 
seconds of data were recorded. This was continued until the transmitted 
power had been reduced by a total of 70 dB. The effect of reducing trans- 
mitted power is to decrease SNR although it may not be in the same steps of 
10 dB because of the action of the rocket receiver AGC. Since the rocket 
was stationary, the output of the rocket receiver only consisted of the 
reference frequency, harmonics of the reference frequency, and noise. By 
using the frequency analysis procedure described in Chapter 5, the measured 
reference signal frequency from the receiver could be compared with that 
measured from the corresponding reference signal channel. The difference 
between these two measured frequencies is an indication of possible error in 
Faraday rate that would result from various SNR levels. 

Figure 6.1 shows the results of this test for system 1 (2.225 MHz) as 
far as the 40 dB level of power reduction. As is indicated, the peak 
magnitude of the random error increased as the transmitted power was reduced. 
This fact is also represented numerically, in Figure 6.1, by the standard 
deviation (a) for the error in the frequency which is calculated for each 
one-second interval. The standard deviation typically increases two to three 
times for each 10 dB reduction in transmitted power. 

During the test just described the ground transmitter was located 
approximately 1 km from the rocket. However, the experiment is primarily 
designed to operate with the rocket at an altitude of the order of 100 km 
(the range is typically 50 to 150 km) . The additional free space propagation 
path loss between 1 km and 100 km is approximately 26 dB. Therefore, it would 
be expected that if the rocket could be placed at a 100 km distance from 
the transmitter and the SNR system test repeated with full power, the standard 
deviation of the error in the frequency measurement would be about 9 mHz; 
that is, from Figure 6.1, between the 20 dB standard deviation level of 
6.0 mHz and the 30 dB level of 10.5 mHz. Further as the rocket moved to 
greater distances it would also be expected that the magnitude of error would 
continue to increase. This type of behavior is indicated in Figure 6.2 which 
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is actual data from the flight of Nike Apache 14.534, system 1. In this 
figure the frequency of reference signal measured from the reference channel 
(i.e., direct line to the TM station) and from the receiver channe 1 
(i.e., through the transmitter/i Jceiver/telemetry system) is shown as a 
function of time after launch. As this figure indicates the frtouency of the 
reference channel remains well-behaved during the flight while the frequency 
measured from the receiver channel increasingly becomes more erratic. This 
is interpreted to be a direct result of decreasing SNR. 

The standard deviations for the measurement of receiver frequency given 

in Figure 6.1 cannot be directly compared with uncertainty in electron density. 

This is because electron density is only indirectly related to the frequency 

measurement through Faraday rotation and equation (5.3). However, since two 

receiver signal frequencies are required in equation (5.3), the standard 

deviations for frequency measurement error can be used in this equation to 

indicate an uncertainty for Faraday rate. By assuming, in equation (5.3), 

that the errors in both measurements of receiver frequency are statistically 

independent, the corresponding error in Faraday rate can be estimated. For 

a standard deviation of 9 mHz, the expected error at 100 km, the uncertainty 

in Faraday rate is 4.5 mHz. By calculation an uncertainty in Faraday rate 

of 5.7 mHz at the 2.225 MHz transmitting frequency corresponds to an 

electron-density uncertainty of 100 cm"^. Therefore, when this uncertainty 

has been combined with other possible errors, it should be expected that the 

electron-density error will be about 100 cm'"*. Figures 5.6 and 5.7 show 

that the measured electron-density values for Nike Apaches 14.533 and 14.534, 

.3 

respectively, both have an uncertainty on the order of 100 cm for the final 
calibrated probe current in the lower portions of the profile (i.e., near 
100 km) . By increasing SNR it is expected that this uncertainty would be 
decreased. 

The SNR system test also provides detailed information about the 
receiver output spectrum as a function of input signal level. In general it 
is found that as the transmitted signal level is reduced, the output signal 
decreases very slowly and the noise level increases rapidly. Table 6.1 shows 
this result numerically by giving the relative amplitude of the reference 
signal as well as the peak noise level observed within a 30 Hz bandwidth 
centered at the reference frequency. The noise level at frequencies near 
that of the reference signal is an indication of the level of noise present 


Table 6. 1 Peak receiver output signal and noise levels 
versus relative input signal level. 


Relative Input 
Signal Power (dB) 

Approximate Reference 
Signal Peak Amplitude 

Peak Noise Amplitude 
Within 30 Hz B.W. 

(Peak Signal] 
Peak Noise 

0 

76 

0.16 

54 

-10 

75 

0.49 

44 

-20 

72 

0.96 

38 

-30 

66 

2.5 

28 

-40 

50 

4.6 

21 

-50 

20 

7.8 j 

8 


Note: 


All amplitudes are relative. 
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at the reference signal frequency. This "in-band" noise is directly 
responsible for frequency measurement error. The information in this table 
was obtained by using the FFT spectral analysis procedure of Section 5.1. 

Table 6.1 shows that the signal spectral amplitude decreased very slowly as 
the input signal power level was reduced in 10 dB steps. This is the result 
of AGC action in the receiver. On the other hand, the peak noise eunplitude 
approximately doubled for each 10 dB reduction of input signal power. The 
ratio of the peak amplitudes of signal and noise is also calculated in decibels 
for each of the input power levels. For the 0 dB relative input level, the 
output signal was approximately 500 times stronger than the background spectral 
noise. This factor decreased as the input signal power was reduced. 

6.2.3 Noise in the system. The second system performance check 
undertaken was designed to monitor the effects of noise at several locations 
within the experimental system. The procedure was to use portable equipment 
to monitor the effect of noise on the instantaneous frequency of the 
reference signal. The equipment for this purpose is shown in Figure 6.3 where 
it is seen that the ief;rence signal is first band-pass filtered (BPF) then 
FM detected followed by a low-pass filter (LPF) and finally spectrum analyzed. 
The concept that is used by this approach is that any noise within the 
bandwidth of the BPF causes the instantaneous frequency to vary which is 
demodulated by the FM detector. Spectral analysis of this detected signal 
allows for both the identification and classification of any noise present; 
that is, determine the presence of any noise and whether it is random or 
periodic. Since the digital data processing provides very narrow bandwidth 
filtering, only noise within several Hz of dc of the FM detector output 
are of concern. 

For this test the experimental system was also operated in a pre- launch 
configuration and measurements were made at the transmitter van, launch pad 
blockhouse, and at the telemetry station. The van measurements were of the 
reference signal prior to amplification of the two CW carriers. The block- 
house measurements were made from the rocket payload control panel which 
monitored the receiver output. In the telemetry station, measurements were 
taken of the signal that was digitally sampled for recording purposes. 

Results for the test can be divided into those related to system 1 and 
system 2. For system 1 no noise was observed at the van but white noise 
was observed at both the blockhouse and the telemetry station. The magnitude 
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of the white noise was inversely related to the transmitted signal amplitude. 
For system 2 white noise was also observed at the blockhouse and telemetry 
station, but not in the van. In addition, periodic noise was also observed 
at 40, 60 and 80 Hz at all three locations. The source of noise for the 40 
and 80 Hz components was identified as coming from the equipment that 
generates the 999.75 and 1000.25 kHz carrier signals in system 2; it was 
observed that the periodic noise frequencies could be moved several hertz by 
slight adjustments to the carrier signal frequencies. It was also discovered 
from this test that the system 2 local oscillator, f , shown in Figure 4 2 
did not have the best type of crystal oven. Instead of using a linear 
(proportional) oven, a switching type was being used. As detailed in Section 
4.1 switching ovens cause transient behavior in the reference frequency; this 
problem has subsequently been corrected. 

Several conclusions can be made from this system test. For system 1, 
the source of the noise could not be clearly identified; from the previous 
SNR system test it is clear that the source was located before the receiver 
output but after the oscillator generating equipment in the van. This 
conclusion also holds for system 2 but further conclusions are possible. 
Theoretically, a constant amplitude sinusoid that frequency modulates a 
carrier generates a spectrum that places sidebands at multiples of the 
sinusoid frequency around the carrier frequency. The FM modulation from the 
40, 60 and 80 Hz frequencies should not interfere with the frequency measure- 
ment of the carrier since the sidebands would be outside the narrow signal 
processing bandwidth. However, since it was demonstrated that the 40 and 80 
Hz components could be shifted significantly in frequency, it is possible 
that they could, under some circumstances, cause interference. 

6.2.4 Conclusions. As a result of the system tests as well as experience 
from the data analysis of Nike Apaches 14.533 and 14.534 certain sections of 
the system can be eliminated as major sources of noise under normal operating 
conditions. To begin with, the telemetry system does not .' -pear to introduce 
any noise. Th..' can be seen by observing the magnetometer rocket roll signal 
which uses the same or identical telemetry and data recording equipment as is 
used by the receiver channel. The rocket roll rate decreases verv *owly once 
the rocket is free of the atmosphere. In the measurement of the rocket roll 
rate, using the spectral analysis techniques cf Section 5.1, this character- 
istic was observed; it is shown in Figure 6.4 for Nike Apache 14.534. This 
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Figure 6.4 Rocket roll rate versus seconds after launch for Nike Apache 14.534. 
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plot is significant because it does not show any effects attributable to 
noise. Therefore, it is concluded that the noise irtroduced in the telemetry 
system, the A-to-D conversion, and the data recording equipment is negligible. 

Some parts of the ground transmitting equipment can also be eliminated 
as noise sources. This is because it would be expected that any noise 
generated in the oscillators, mixers, or amplifiers would be most apparent 
at full power. However, the results of the SNR test of system 1 with full 
power indicate that no noise of significant amplitude exists. No similar 
conclusion can be directly made about system 2 from the tests, but it does 
demonstrate that equipment of this complexity is adequate for experiment 
requirements, when operating correctly. Further, the analysis of the 
reference t gnai channel also indicates no problem with either system with 
the exception of the frequency jump problem in system 2 due to the 
switching -type crystal oven. 

The combined results and conclusions of the two system tests above 
indicate that the major source (or sources) of noise involve the receiver 
and atmospheric noise. It is very likely that ?.t these frequencies of trans- 
mission and for moderate signal levels the dormant noise sources wou^d be 
atmospheric "sky" noise, however. Figure 4.4 indicates that the receiver can 
coi. tribute significant amounts of noise ' jr low signal levels. Further, the 
SNP test has shown that at sufficiently large signal levels the experimental 
system is able to operate with low random error. 

63 Measurement of Frequency by Spectral Analysis 

Analysis of system performance is extended in this section to include 
the computer frequency measurement process of digital data. Several factors 
have been identified that presently limit the overall accuracy of this 
procedure. These limitations, which are both theoretical and practical in 
nature, deal with the definition of average frequency, the problem of 
frequency measurement overshoot, the effect of Gaussian windowing, and 
inaccurate rocket roll rate measurements. 

6.3.1 Step changes in frequency . Section 5.1 presented a Fourier 
spectral analysis technique for measuring the required frequencies in 
equation (5.3) to obtain Faraday rotation rate. This procedure requires one 
second of data for each frequency measurement; taken. In general, the 
signals that are measured for equation (5.3) not have a constant frequency 
during the measurement interval and, therefore, some ambiguity exists in 
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interpreting the resulting measured values. It would be hoped that the 
values returned by the measurement process would correspond to the frequency 
time average for that interval. However, it has been established that this 
is generally not the case. 

As an example, the theoretical spectrum of a one-second signal was 
calculated and the frequency of peak amplitude determined. The signal was a 
sinusoid at 490.50 Hz from t=0tot=0.5s and then suddenly changed in 
frequency to 490.51 Hz for the remainder of the interval: t » 0.5 s to 

t = 1.0 s. The peak spectral amplitude was found to be at 490.512 Hz. On 
the other hand, the time-averaged frequency for this signal was only 490.505 
Hz. This example shows that the frequency cf a spectral peak need not agree 
with the time-averaged frequency for a signal. 

Theoretically, it is expected that the spectral peak frequency and the 
time-averaged frequency should agree only for a constant frequency sinusoid 
or a sinusoid whose frequency varies in a purely random fashion [ Taut and 
Schilling, 1971]. For an actual electron-density profile the above results 
imply that reliable data can only be obtained in those regions with constant 
electron density or those with complex structure. In practice, regions of 
slowly varying electron density also have been found to give acceptable 
results. However, for other regions of the profile this ambiguity between 
measured frequency and actual time-averaged frequency can be a cause of 
error. 

6.3.2 Varying frequency. An additional problem has been observed in 
measuring frequencies which vary during the measurement interval. It has 
been found that when steep gradients of electron density commence within a 
measurement window, the measured frequency tends to significantly overshoot 
or undershoot the actual frequency. Figure 6.5 shows the results of an 
example that will help clarify this problem. Artificial data was generated 
from a sinusoid which undergoes an interval of linearly increasing frequency 
as seen in the figure. The slope of this linear increase in frequency was 
such that a 30 mHz change occurred in a half second. 

As the window was slid across the region of changing frequency the 
measured frequency significantly overshot the actual value which is also 
indicated in Figure 6.5. The maximum error between measured and actual 
frequencies was approximately 66 mHz even though the total actual frequency 
change was 50 mHz. The corresponding time-avtr? 2 ed value agreed much more 
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closely as also seen in the figure. 

It has been similarly established that a negative frequency slope would 
v cause the measured frequency to undershoot the actual value. The magnitude 

of this frequency change over a half second is typical of what would be 
expected to occur at ledges in electron-density profiles. However, the 
example above may be artificial in that no attempt was made to prevent 
discontinuities from occurring at the boundaries of the changing frequency 
interval. Discontinuities in data cause spectral energy to spread over a 
wide frequency range which can introduce additional err >r into the measure- 
ments. In general, actual data will not have discontinuities such as these. 

Even so, evidence of overshoot and undershoot effect can be found in the 
experimental results of Figures S.6 and 5.7. 

6.3.3 Windowing. It was also indicated in Chapter 5 that Gaussian 
windowing was needed to reduce the effects of spectral leakage, however, this 
process creates a new problem. This problem results from the fact that the 
amplitude of the Gaussian window function approaches zero near the ends of 
the measurement window. During data processing that data which is aligned 
with the ends of the window interval are significantly attenuated due to 
multiplication with the Gaussian function. The result of this is that the 
frequency spectrum calculated for each interval is primarily determined from 
that da-a in the middle of the interval. Thus the use of windowing functions, 
such as the Gaussian, introduces the additional complication that some data are 
lost or greatly attenuated by the windowing process. This is shown in 
Figure 6.6(a) where the regions of lost data have been indicated. This fact 
prevents the measured frequencies from truly reflecting the data of the 
entire window interval. 

It is of value at this point, for analysis purposes, to present one possi- 
ble correction to this problem of lost data. This is done, at the expense of 
additional computer processing, by inserting new window intervals between the 
original. Figure 6.6(b) shows these new windows with dotted lines which now 
overlap the old windows. As is indicated in the figure the new windows would 
primarily use that data which was lost to the original non-overlapped windows. 

To prevent the number of subsequent computer processing steps from also increas- 
ing, after Faraday rotation calculation, the new Hita .an ee averaged into the 
olu using a process known as Hanning filtering. The Hanning filter is essential- 
ly a method of calculating running averages. For this appli ation the algorithm 
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Figure h.6 l.xample of window overlapping to provide better use of available data. 
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appears as 


(fj) - | ♦ 1 f, * 1 f. 

J FILTERED J J J 


where F. is the j'th Faraday rotation value measured from the j’th original 

3 

non-overlapping window and F. and F. represent Faraday rotation values from 

j - J* 

the new overlapped windows just before and after the jth original window, 
respectively. Using the above algorithm the number of filtered Faraday 
rotation values is equal to the number generated for the original non-over lapped 
windows. 

The use of this alternate procedure fo" calculating Faraday rotation has 
two basic advantages. Firstly, as discussed above the use of overlapping 
windows enables all data to be considered in determining Faraday rotation. 
Secondly, Harris [1978] has indicated that measurement variance due to 
random noise can be reduced by the use of Gaussian windows with 50 percent 
overlap when the overlapped data is averaged. This latter point is signifi- 
cant because as indicated in Section 6.1 random atmospheric and receiver noise 
is the major source of error introduced by the experimental system. There- 
fore, data collected using this method shculd more closely represent the 
actual average Faraday rotation desired as well as being less sensitive to 
the effects of random noise. 

Figure 6.7 provides a means of comparison between the use of non-overlapped 
and the overlapped Hanning filtered measurement procedure. This figure 
shows both methods, using Faraday rotation data from system 1 for the flight 
of Nike Apache 14.533. The results seem to substantiate the expected 
performance: the data derived from the new method ♦ aids to have smoother 
variation between data points when compared with the oid method. This 
behavior is particularly noticeable beyond 90 SAL where the data becomes 
increasingly noisy (random), a predicted response from the previous section 
due to degrading SNR which is further substantiated by this figure. As is 
shown the new method data does not have as large peak variations when Cleared 
to the old. Before 90 SAL the data from both methods agree much more closely 
demonstrating the point that the new method does not simply smooth the old. 

It is important to remember tnal the new method uses three measurement, values 
as opposed to one to generate each data point on the graph. 


ife 
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6.3.4 Rocket roll rate. The -final problem that has been associated 
with the software data processing is in the measurement of the rocket roll 
rate. This measurement is particularly important because it is multiplied 
by 4 in equation (5.3) to determine Faraday rate. Although Faraday rate 
should be zero below the ionosphere, it has been observed for both Nike Apaches 
14.533 and 14.534 (Figures 5.3 and 5.4, respectively) that a periodic error 

is present (prior to approximately 67 SAL). The error oscillates with a 
maximum deviation which is typically about 6 degrees per second and is 
independent of the system measurement frequency (2 or 3 MHz). It is clear 
that the error is not associated with any ionospheric effect and since the 
errors from both systems 1 and 2 are well correlated, the error must be 
caused by a common factor for both systems. It car be further seen in the 
figures that the oscillating error has a period which slowly increases until 
ionospheric penetration where any possible continuation of the error becomes 
obscured by Faraday rotation. This increasing period behavior appears to t> 
related to rocket precessional effects. In general, while the rocket is still 
in the atmosphere, it has a precessional motion which has a period that 
slowly increases. Upon leaving the atmosphere (at about 80 km), the pre- 
cessional period ceases to increase and stabilizes to a nearly constant value. 
The final values vary from flight to flight apparently dependent on the 
behavior of the rocket at the moment of exit from the atmosphere. Precessional 
motion can be observed on the magnetometer signal, which is used to measure 
rocket roll rate, as an amplitude modulation of a sinusoidal signal; it is 
also possible that precessional motion effects the frequency of the sinusoid 
as well. In either case since precession is present throughout the flight 
and since precessional caused error has been observed prior to ionospheric 
p.netration, it is probable that error is being introduced into electron- 
density measurements due to precessional motion. 

6.3.5 Conclusion. Several sources of error in the frequency measure- 
ment relate! to software have been identified in this section. Some, such as 
the definition of average frequency and the problem of frequency overshoot, 
are particularly hard to overcome because they are an inherent part of this 
measurement process. On the other haau the problem of lost data due to 
windowing can be overcome with the overlapping procedure; this method is also 
more tolerant to noise. The error due to precessional motion may also be 
correctable. It would appear that improving system SNR and correcting for 
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precessional motion effects would be particularly beneficial since SNR is 
critical and since rocket roll rate errors are magnified in equation (5.3). 
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7. CONCLUSIONS AND RECOMMENDATIONS FOR FURTHER WORK 

The original reason for modifying the radio propagation experiment was 
to provide an experimental system capable of nighttime operation; the original 
experiment had not been satisfactory at night. The result of the work that 
followed was an experiment capable of providing data of sufficient accuracy 
during the nighttime, with the added advantage of reducing the complexity of 
the system. This was verified for the flights of Nike Apaches 14.533 and 
14.534. However, it was found on those flights that the original experiment 
was also capable of providing useful nighttime data. No conclusive explana- 
tion has been found to explain the original experiment's recent success as 
opposed to its previous failures. It is possible that measurements on 
previous flights were made incorrectly or that ionospheric conditions were 
better fo. the two recent flights. 

The accuracy of the propagation experiment appears to be adequate for its 
present application: the calibration of the Langmuir probe experiment for 
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electron-density values of 100 cm and greater. There are possibilities for 
further improvement. Based on system performance tests, significant improve- 
ments can be expected by improving SNR at the output and by correcting for 
precessional motion. Further improvement can be expected by dealing with 
the other error sources described in Chapter 6. 

Several approaches exist for improving SNR. The first method would be 
to increase transmitted signal power. A minimum increase in signal power of 
10 dB should be attempted; this will reduce the standard deviation of frequency 
measurement error by at least one half, as discussed in Section 6.1. Other 
less costly methods can be used to minimize noise and its effect. For 
example, the window overlap method described in Section 6.2 appears 
particularly attractive for noise suppression and complete data utilization. 
Also noise generated by the payload receiver should be considered. It should 
be determined whether the amount of noise presently generated by the receiver 
is significant when compared with typical atmospheric noise levels. If it 
is, then the receiver can possibly be tuned for minimum noise figure, a 
procedure not presently followed. 

Various other modifications can be considered to improve the accuracy 
and reliability of the experiment. For example, at present no real-time 
pre-launch system check exists that can verify the operational readiness of 
the experiment. The system check would involve confirming adequate signal 
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strength, frequency stability, spectral purity, and the proper operation 
of data recording equipment. Semi-portable equipment to monitor the first 
three of the above items would be particularly helpful for trouble-shooting. 

Certain modifications can be directed towards improving spectral leakage 
interference that occurs during Fourier spectral analysis. As described in 
Section 5.1 spectral leakage occurs whenever a signal is not periodic within 
the window interval (of one second), projecting energy over all spectral 
frequencies. Use of window functions, such as the Gaussian, are designed to 
reduce this problem. However, the situation can be further improved by 
reducing the amplitude of any signal which is potentially not periodic within 
the window interval; this would include any noise or sinusoid that is not 
exactly at an integer frequency (1, 2, 3, ... Hz). A simple modification for 
this purpose is to band-pass filter the receiver output signal prior to FFT 
analysis. This could be accomplished by analog or digital methods. 

Another technique directed towards the elimination of unwanted inter- 
modulation signals is the use of a square law detector instead of the 
envelope detector of the present payload receiver. The desired mixing 
operation of the detector input signal would be undertaken by a squaring opera- 
tion. The higher order terms present in the output of the envelope detector 
generate unwanted components that have the potential of causing spectral 
leakage interference. A square law detector would not have the higher order 
components present. Neither of these suggested modifications would requite 
major system alterations. Relatively simple checks could be devised to verify 
the usefulness of these changes before actually carrying out the modification. 

Another problem that should be considered is that of impulse noise, 
atmospheric noise due to static discharges from lightning can generate 
impulsive radio interference that can propagate many hundreds of miles at HF. 
This impulse noise is often observed in the received signal, particularly at 
low signal levels. The exact effect of impulse noise on FFT frequency measure- 
ments has not yet been determined, but it should be expected that an impulse 
will generate a ”hite noise power spectrum, thus increasing the error. Some 
kind of limiter action in the receiver may be able to reduce the effect. 

It is hoped that some of the suggested modifications above can be imple- 
mented for future nighttime rocket flights. It is expected that improved 
accuracy and more reliable operation can be obtained whi’e still maintaining 
overall experimental simplicity. 
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APPENDIX I. PROGRAM FOR READING DIGITAL TAPES 


FRIES. 

SI3NON' 1 225»155:<). 

mmm - — (code woe ? core here) 

61 Li. , cLEC , PS6399 • 

PRINT, RJE=EE/CC/EJ. 

SETTL , 30 . 

FTN,A,L=0. 

LABEL , TAPE 1 , NT , LB=KU , P0=R , VSN=WI54 3F-A047 , FsL . 

ATTACH , UOILIB/UNsLIBRARY . 

LD3ET(LIB=SYMLIB/U0ILIB) . 

LGO. 

SKIPTO.X. 

EXIT. 

ENDSKIP.X 

DAYFILE. 

?'~ opa >! (input,output,tapei , tape2=output , tape3-input ) 

* THIS PROGRAM IS INTENDED TO DISPLAY THE HEADER RECORDS AND DATE 

* CODES AND/OR CHECK THE TIME BETWEEN RECORDS ON THE DATA TAPE 

* AGAINST THE EXPECTED VALUE FOR NASA DATA TAPES RECORDED AT WALLOPS 
» ISLAND VIRGINIA. THE USER SUPPLIES THE BEGINNING RECORD NUMBER 

* AND THE NUMBER OF RECORDS TO BE READ INDICATING IF ALL TIME RECORDS 

* IN THAT BLOCK WILL BE PRINTED OR JUST THOSE WITH TIME LENGTHS NOT 

* EXPECTED. MORE THAN ONE DATA CARD CAN BE USED. DATA CARDS MUST 
» HAVE THE FOLLOWING FORMAT: 


FIRST DATA CARD 

COLUMN 1-4 
OTHER DATA CARDS 
COLUMN 1-4 
COLUMN 5 
COLUMN 5-9 
COLUMN 10 
COLUMN 11 


EXPECTED TIME BETWEEN RECORDS, MILLISECONDS* 1 0 

FIRST RECORD IN BLOCK, RIGHT JUSTIFIED 
COMMA 

NUMBER OF CONSECUTIVE RECORDS TO CHECK 
COMMA 

1 IF ALL TIME CODES TO BE PRINTED 
0 IF ONLY DESCREPENCIES PRINTED 


» A 789 CARD MUST FOLLOW THE DATA DECK. THE PROGRAM WILL PRINT HEADERS 

* FOR RECORDS OF LENGTH 10, 45, 1005, OR 2008 WORDS. EACH WORD IS 16 

* BITS LONG. 

* WRITTEN BY KEITH FRIES JULY 5, 1978 

» AERONOMY LABORATORY - UNIVERSITY OF ILLINOIS 
IMPLICIT INTEGER (A-Z) 

REAL DIFF.DFRSEC,F1 

DIMENSION ARRAY ( 2008 ) , PARRA Y ( 400 ) 

DATA (NBLOCK=0) 

* THE FOLLOWING LINE LOCATES THE 'FILE ENVIRONMENT TABLE' ADDRESS WHICH 

* IS USED LATER IN THE PROGRAM TO BACKSPACE TAPE1 . 

NN=FETADR( 1 ) 

* SET FLAG FOR PRINTING HEADE: WITH FIRST GROUP OF RECORDS 

HEAD=0 

* WRITE "NASA" BANNER AT TOP OF OUTPUT. 



WRIT£(2,5) 


. FORMAT! "0"/"0",6X 

»"N 


N AAAA 


SSSS AAAA" , 1 3X , "AERONOMY" 

J/" ",6X,"NN 

N 

A 

A 

S S 

A 

A" 

J/" ",6X,"N N 

N 

A 

A 

s 

A 

A", 10X, "LABORATORY" 

.1/" ",6X,"N 

N N 

AAAAAA 

SSSS 

AAAAAA" 

J/" ",6X,"N 

N N 

A 

A 

s 

A 

A", 14X, "ROCKET" 

J/" ",6X,"N 

NN 

A 

A 

S S 

A 

A" 

J/" ",6X,"N 

N 

A 

A 

ssss 

A 

A", 13X, "PROGRAM" 

J/"0",11X,"-~ 

CDC 

CYBER 

VERSION — 1 

"/"0 

l"/" 0 ") 


» READ CARD INFORMATION. 

» INPUT STANDARD TIME BETWEEN RECORDS 
READ (3, 6) STDSTP 
6 FORMAT (14) 

• INPUT DESIRED RECORD INFORMATION 

10 CONTINUE 

RE AD (3, 20) RECNUM , NUM , ALLPRT 

20 FORMAT! 14 , IX, 14, IX, II) 

IF(E0F(3) ) 200,21 

•SET FLAG INDICATING FIRST RECORD IN GROUP 

21 START=0 

• PRINT NASA WITH FIRST DATA GROUP, ELSE START NEW PAGE 

IF(HEAD.SQ.O) GOTO 300 
WRITE(2,310) 

310 FORMAT! "1") 

300 CONTINUE 
HEAD=1 

WRITE! 2, 22) NUM, RECNUM 

22 FORMAT! "0"/"0", "SEARCH ",I4," RECORDS BEGINNING WITH RECORD " 
J ,14 , " .") 

IF( ALLPRT. EQ. 1 ) GOTO 24 
WRITE(2,23) 

23 FORMAT!" "."ONLY RECORDS WHICH HAVE NON-STANDARD TIME ", 
J"BETWEEN RECORDS WILL BE PRINTED") 

24 CONTINUE 
WRIT£(2 ,26) 

26 FORMAT! "0"/" ") 

• SKIP THE PROPER NUMBER OF RECORDS TO GET TO THE RECORDS DESIRED. 

25 ISKIP=RECNUM-NBLOCK- 1 
IF(ISKIP) 30,70,50 

30 ISKIP=-ISKIP 
DO 40 1=1 , ISKIP 
40 CALL BACKUP (NN.NBLOCK) 

GOTO 70 

50 DO 60 1=1, ISKIP 
60 CALL FORWRD! 1 .NBLOCK) 

• FILL ARRAY WITH DATA. 

70 DO 190 1=1, NUM 

CALl. TPGET ( 1 , ARRAY , 1 , NWORDS , NBLOCK ) 

IF(NW0RDS . EQ . 2008 ) GOTO 90 
IF ( NWORDS. EQ. 1005) GOTO 130 
IF(NWORDS.LE.45J GOTO 160 
WRITE(2,80) NBLOCK, NWORDS 

30 FORMAT "RECORD ",I4," NON-STANDARD LENGTH, ",I4, 

J" WORDS. "/"0") 
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GOTO 190 
90 CONTINUE 

CALL CALTIM( ARRAY ( 2006 ) , ARRAY ( 2007 ) , ARRAY ( 2008 ) , DAY , HOUR , MIN , SEC , 
JFRCSEC) 

DUDAY/100 

D2= (DAY-D 1*100)/ 10 

D3= (DAY-D 1*100-D2*10) 

HUH0UR/10 

I12=HOUR-H1*10 

MUMIN/10 

M2=MIN-M1* 10 

S1=SEC/10 

S2rSEC-S1*10 

F1=FRCSEC*.0001 

• CONVERT TIME TO SECONDS 

• KEEP IN INTEGER VARIABLES TO AVOID ONVERSION ERROR 

DECSEC=864000000*DAYt-36000000*HOUR>600000*MIN+10000*SRG+FRCSEC 

IF(START.EQ.O) GOTO 126 

IDIFF=DECSEC-OLDDEC 

• CHECK TIME INTERVAL, IF CORRECT DO NOT PRINT 

IF( IDIFF . EQ . STDSTP ) GOTO 126 

• BREAK INTO SINGLE DIGIT VARIABLES IN ORDER TO DISPLAY 

• LEADING ZEROS 

• CONVERT TO REAL TO SIMPLIFY PRINT ROUTINE, CONVERSION ERROR 

• MILL BE ROUNDED VIA PRINT 

DIFF=IDIFF* .0001 

DDAY=DIFF/86400 

DHOUR= ( DIFF-DDAY*86400) /3600 

DH 1 =DH0UR/ 10 

DH2=DH0UR-DH 1*10 

DMIN=(DIFF-DDAY*86400-DHOUR»3600)/60 

DM1=DMIN/10 

DM2=DMIN-DM1*10 

DSEC=DIFF-DDAY*86400-DHOUR*3600-DMIN*60 
DSUDSEC/10 
DS2=DSEC-DS 1*10 
DFRSEC=DIFF-INT(DIFF) 

TEMP=NBL0CK-1 

WRITE (2,1 24 )NBLOCK, TEMP, DDAY.DH1 ,DH2,DM1 ,DM2,DS1 ,DS2,DFRSEC 
124 FORMAT("+",60X,"TIME BETWEEN RECORDS ",I4," AND ",I4, 

J": DAYS " , 1 3 , " , " ,21 1 , " : ” ,21 1 , " : " ,21 1 ,F5 . 4 ) 

GOTO 127 

126 CONTINUE 

IF( ALLPRT.NE. 1 ) GOTO 129 

127 CONTINUE 

WRITE(2 , 1 10)NBLOCK,D1 ,D2,D3,H1 ,H2,M1 ,M2, S 1 ,S2, FI 
110 FORMATO’-*-", "RECORD ",I4," TIME: DAY ", 

J3I1 " , 21 1 , " : " ,21 1 , " : ,21 1 ,F5 . 4 , " U.T.'V'O") 

129 CONTINUE 
START= 1 
OLDDECrDECSEC 
GOTO 190 

130 WRITE(2, 140) NBLOCK 

140 FORMAT("+’\ "RECORD ",I4," CONTAINS 1005 WORDS. "/"O") 


- -*■*”«* "»•>* 


GOTO 190 

160 WRIT£(2, 170) NBLOCK , NWORDS 95 

170 FORMAT("+", "RECORD ",I4," CONTAINS ",I2," WORDS.") 

• THE FOLLOWING THREE DO-LOOPS ARE NECESSARY TO CONVERT TEXTUAL DATA ? 

• FROM IBM EBCDIC REPRESENTATION TO CYBER DISPLAY REPRESENTATION. > 

DO 173 N=1, NWORDS 
M=2*N-1 

173 CALL GBYTES(ARRAY(N),PARRAY(M),44,8,0,2) 

M=2*NWORDS 
00 176 N=1.M 

176 PARRAY(N)=CONETD(PARRAY(N) ) 

DO 178 N=1 ,M 

IF ( ( PARRAY(N) . EQ. 00) .OR . (PARRAY(N) .EQ. 42) .OR . (PARRAY(N) . EQ.50) .OR . 

J(PARRAY(N) .EQ.63) ) PARRAY(N)=45 
173 CONTINUE 

WRITE(2, 130) (PARRAY(N) ,N=1 ,M) 

1 30 FORMATC" " ,90R1 ) 

WRITE(2 , 182) 

132 FORMA'K "0" ) 

190 CONTINUE 
GOTO 10 
200 STOP 
END 

SUBROUTINE CALTIM( T2006 , T2007 , T2008 , DAY , HOUR , MIN , SEC , FRCSEC) 

IMPLICIT INTEGER (A-Z) 

DIMENSION 1(16) 

••SUBROUTINE CALTIM (CALCULATE TIME) CALCULATES TIMES 
••BY DECODING THE NASA DIGITAL TIME DATA IN THE LAST 
«• THREE WORDS OF A TAPE RECORD. THIS IS UNIVERSAL 
••COORDINATED TIME IN DAYS, HOURS, MINUTES, SECONDS, 

• AND DECIMAL FRACTIONAL SECONDS. 

• CONVERT T2006 INTO FRACTIONAL SECONDS. 

CALL BINARY(T2006,I) 

FRCSEC=(8*I( 16)+4*I( 15)+2*I( 14) +1(13) )*1000 
1+(8*I(12)+4*I(11 )+2*I( 10)+I(9))*100 
2+(8*I(8)+4*I(7 )+2*I(6)+I(5) )*10 
3+(8*I(4)+4*I(3)+2*I(2)+I( 1 ) ) 

• CONVERT T2007 INTO MINUTES AND SECONDS. 

CALL BINARY (T2007, I) 

MIN=(4*I( 15)+2*I( 14)+I(13))*10 
1+( 8* 1(12) +4*1(11 )+2*I( 10) +1(9)) 

3EC=(4*I(7 ) +2*1 (6) +1(5) )*10 
1+(8*I(4)+4*I(3)+2*I(2)+I( 1)) 

» CONVERT T2008 INTO DAYS AND HOURS 
CALL BINARY ( T2008 , I ) 

DAY=(2*I( 16)+I( 15) ) * 1 00 
1+(8*I( 14) +4*1 C 1 3)+2*I( 12)+I( 11 ) )*10 
2+(3*I ( 10) +4*1 ( 9 ) +2*i ( 8 ) +1 ( 7 ) ) 

HOUR= (2*1(6) +1(5) )*10 
1+( 3*1(4) +4*1(3) *-2*1 (2) +1(1)) 

RETURN 

END 
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SUBROUTINE BINARY ( IINTGR , I ) 

INTEGER 1(16) 

• SUBROUTINE BINARY GENERATES THE BINARY REPRESENTATION 

• OF THE NUMBER INTGR, FOR USE IN SUBROUTINE CALTIM. 

INTGR=I INTGR 
DO 1 N=1,16 

1 I(N)=0 

DO 2 N=1 , 16 

J=17-N 

L=2«(J-1) 

IF( INTGR. GE.L) I(J)=1 

IF( INTGR. GE.L) INTGR=INTGR-L 

2 CONTINUE 

• THE ARRAY 1(1-16) IS THE BINARY REPRESENTATION OF 

• THE NUMBER INTGR, WITH 1(16) BEING THE MOST 5IGNI- 

• FICANT BIT,2»*15, AND 1(1) THE LEAST SIGNIFICANT 

• BIT,2**0. 

RETURN 

END 

SUBROUTINE TPGET (U, ARRAY, M,N, NBLOCK ) 

• THIS SUBROUTINE IS USED TO TRANSFER ONE RECORD OF TAPE DATA FROM 

• UNIT »U' TO 'ARRAY', BEGINNING WITH ARRAY ELEMENT M. THE NUMBER 

• OF lo-BIT WORDS TRANSFERRED ON EACH CALL IS N (TYPICALLY 10, 45, 

• 1005, OR 2008). THE NUMBER OF THE RECORD JUST READ IS NBLOCK. 

IMPLICIT INTEGER(A-Z) 

DIMENSION BUFFER (540), ARRAY ( 1 ) 

BUFFER IN (U, 1 > (BUFFER ( 1 ) .BUFFER (540) ) 

NBL0CK=NBL0CK+1 
IF(0NIT(U) ) 10,20,40 
10 CALL LENGTHX(U,I, J) 

• I IS THE NUMBER OF 60-BIT WORDS READ. 

• J IS THE NUMBER OF BITS IN THE LAST 60-BIT WORD THAT WERE NOT USED. 

N=(60 # I-J)/16 

CALL GBYTES ( BUFFER , ARR A Y ( M) , 0 , 1 6 , 0 , N ) 

RETURN 

20 NBLOCK=NBLOCK- 1 
WRITE (2, 30) NBLOCK 

30 FORMAT (" "/”0","**« END OF FILE ENCOUNTERED. LAST RECORD IS ",I4, 
J". »**") 

STOP 1 

40 WRITE (2, 50) NBLOCK 

50 FORMAT ( " »/"0"," # ** PARITY ERROR DETECTED IN RECORD ",I4,". ***") 
STOP 2 
END 


SUBROUTINE FORWRD(U, NBLOCK) 
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? * THIS SUBROUTINE IS USED TO FORWARD SPACE ONE RECORD ON TAPE UNIT U. 

• THE NUMBER OF THE RECORD SKIPPED IS NBLOCK. 

IMPLICIT INTEGER(A-Z) 

DIMENSION 3UFFERC540) 

BUFFER IN CU, 1 ) ( BUFFER ( 1 ) , BUFFER (540) ) 

NBLOCK=NBLOCK+ 1 
IF(UNIT(U) ) 10,20,40 
10 RETURN 

20 NBL0CK=N8L0CK- 1 
WRITE (2, 30) NBLOCK 

30 FORMAT (" "/"0","*«* END OF FILE ENCOUNTERED. LAST RECORD IS ",I4, 
J". **»") 

STOP 3 

40 WRITE( 2 ,50) NBLOCK 

50 FORMAT (" PARITY ERROR DETECTED IN RECORD ",I4,". •«•") 

RETURN 
END 


INTEGER FUNCTION FETADR(UNIT) 

• THIS INTEGER FUNCTION RETURNS AS ITS VALUE THE ADDRESS OF 
» THE "FILE ENVIRONMENT TABLE" CORRESPONDING TO THE FORTRAN 

• UNII NUMBER GIVEN BY "UNIT", OR ELSE 0, IF THAT UNIT WAS NOT 

• DEFINED IN THE "PROGRAM" STATEMENT IN THE MAIN PROGRAM. FOR 

• INSTANCE, IF THE PROGRAM STATEMENT IS 

» PROGRAM XYZ(TAPE1 ,TAP’C2, INPUT, TAPE5=INPUT) 

• THEN FETADR(I) WILL BE THE MACHINE ADDRESS OF THE "FET" FOR 

• "TAPE1 " , FETADR(2) WILL BE ADDRESS OF THE FET FOR TAPE2 

« AND FETADRC5) WILL BE THE ADDRESS OF THE FET FOR FILE "INPUT", 

• SINCE UNIT 5 IS EQUATED TO "INPUT" IN THE PROGRAM. FETADR(6) , 

• FOR EXAMPLE, WILL BE 0, SINCE TAPE6 DOES NOT APPEAR IN THE 
« PROGRAM STATEMENT. 

» THIS FUNCTION IS INTENDED TO WORK ONLY UNDER CDC CYBER FORTRAN. 

• IF YOU THY TO USE IT IN ANY OTHER CIRCUMSTANCE, YOU DESERVE 

• WHAT YOU GET. ALSO, UNLESS YOU KNOW WHAT A FET IS AND WHAT IT'S 

• FOR, DON'T MUCK WITH THE ROUTINE, BECAUSE OTHER ROUTINES MAY 

• DEPEND ON THIS ONE. 

• CODED FEBRUARY 1973 BY STAN KERR AT THE COMPUTING SERVICES OFFICE 

• OF THE UNIVERSITY OF ILLINOIS AT URBANA, FOR BOB ZIMMERMAN OF 

• ELECTRICAL ENGINEERING. 

IMPLICIT INTEGER(A-Z) 

DIMENSION C0RE(1) 

IFfUNIT .GE. 10) GOTO 50 

LFN = 4LTAPE + SHIP I ( 1R0+UNIT, 30) 

GOTO 60 

50 LFN = 4LTAPE + SHIFT ( 1R0+UNIT/10, 30) + 


• > * 



* SHIFT( IR0+M0D(UNIT,10),24) 

60 B = L0CF(C0RE(0) ) 

FETADR = 0 

DO 80 I = 2,64 

IF(CORE(I-B) .EQ. 0) GOTO 100 
1F( (CORE( I-B) .AND. -777777B1 .EQ. LFN) GOTO 90 
80 CONTINUE 

90 FIT ADR = CORE(I-B) .AND. 777777B 

FETADR = C0RE(FITADR-B+1 ) .AND. 777777B 

100 RETURN 

END 


SUBROUTINE BACKUP ( FETADR, NBLOCK) 

THIS ROUTINE IS INTENDED MAINLY FOR BACKING UP 1 BLOCK ON 
A MAGNETIC TAPE BEING READ BY A FORTRAN PROGRAM. IT IS ASSUMED 
THAT THE TAPE IS ACCESSED BY THE PRC ISAM VIA SOME UNIT NUMBER, 

SAY 5. THE USER OF THIS ROUTINE SHOULD FIRST USE THE INTEGER 
FUNCTION "FETADR" TO ESTABLISH THE ADDRESS OF THE "FILE 
ENVIRONMENT TABLE" FOR THE TAPE, AND SAVE THE VALUE (NEVER MIND 
WHY, JUST DO IT). THEN THIS ROUTINE MAY BE CALLED AT WILL TO 
BACKSPACE THE TAPE BY ONE BLOCK. THE SECOND PARAMETER IS SUPPOSED 
TO BE AN INTEGER VARIABLE WHICH THE USER PRESUMABLY USES TO 
KEEP TRACK OF WHICH BLOCK THE TAPE IS POSITIONED ON; ALL BACKUP 
DOES WITH "BLOCK" IS SUBTRACT ONE FROM IT. 

IF THE TAPE IS UNIT 5, THEN THE PROCESS OF USING BACKUP WOULD BE 
SOMETHING LIKE THIS: 

IADDR = FETADRC5) 


CALL BACKUP (IADDR, NBLOCK) 


DO 25 I = 1,10 

25 CALL BACKUP (IADDR, NBLOCK) 

ETCETERA , ETCETERA , ETCETERA 

THIS ROUTINE IS INTENDED TO RUN ONLY UNDER CYBER FORTRAN AT THE 
UNIVERSITY OF ILLINOIS. IT USES AN EXTERNAL ROUTINE "3KIPB" FROM 
SYSTEM LIBRARY SYMLIB TO DO THE ACTUAL BACKSPACE. 

CODED BY STAN KERR AT THE COMPUTING SERVICES OFFICE OF THE UNIVERSITY 
OF ILLINOIS AT URBANA FOR BOB ZIMMERMAN OF ELECTRICAL ENGINEERING. 

IMPLICIT INTEGER(A-Z) 

DIMENSION CORE(I) 
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BIAS = LOCF(CORE(0) ) 

IN = CORE (FETADR+1 -BIAS) .AND. 777777B 

C0RE(FETADR+2-BIAS) = (CORE(FETADR+2-BIAS) .AND. MASK(42)) + IN 
CORE ( r ETADR+3-3IAS ) = (CORE(FETADR+3-BIAS) .AND. MASK(42)) + IN 

CALL SKIPB(CORE( FETADR-BIAS) ,1,1) 

NBLOCK = NBLOCK-1 

RETURN 

END 


IDENT CONETD 
( ( ENTRY CONETD 

*» CONETD - TABLE ORGANIZATION 

* THE CONVERSION TABLE USED BY CONETD IS SET UP WITH 

* 8 DISPLAY CODE CHARACTERS PER WORD, LEFT-JUSTIFIED, TO 

» 3IMPLIFIY THE ARITHMETIC NECESSARY TO ACCESS A WORD 

* WHICH CONTAINS THE DISPLAY CODE CORRESPONDING TO A GIVEN 

« EBCDIC CHARACTER. 


CONETD EQ 
SA 1 
SX5 
NG 
SX2 
1X2 
PL 
BX2 
AX2 
SA2 
MXO 
BX1 
5X1 
LX1 
BX3 
LX3 
1X1 
3B2 
LX2 
MXO 
BX6 
EQ 


0 

XI 

-1 

XI .CONETD 
256 
XI -X2 
X2, CONETD 
XI 
3 

X2+TABLE 

57 

-XO»X1 
X1 + 1 
1 

XI 

1 

XWX3 

XI 

B2 

54 

-X0*X2 

CONETD 


XI = I 

X6=-1 , DEFAULT VALUE IF I OUT OF RANGE 
RETURN IF I<0 


RETURN IF I>255 
X2=I/8 

X2 = TA3LE WORD WITH CONVERTED VALUE 

XI =MOD( 1,8) 

XI =MOD( I ,8)+1 


X 1 =6* C 1+MODCI.8) ) 

SHIFT CONVERTED VALUE TO LOW 6 BITS 
AND OUT EXTRANEOUS AND RETURN 


TABLE 


DATA 55626460555255650000B,55555545564657500000B 
DATA 333435365555515500008 , 434455556154737 10000B 
DATA 55555555554742770000B.55555555550067700000B 
DATA 55554 155555555530000B, 555555553740550000005 
DATA 555555555555555500008,555561 57725 145660000B 
DATA 67555555555555550000B.55556253475277760000B 
DATA 46505555555555550000B, 555575566365737 10000B 
DATA 55555555555555550000B, 557400607 47054640000B 
DATA 55010203040506070000B, 101 15555555555550000B 
DATA 55121 314151617200000B.21225555555555550000B 



100 


DATA 5576232U252627300000B , 3 1 325555555555550000B 

DATA 55555555555555550000B , 55555555555555550000B 

DATA 720 1 0203040506070000B ,1011 5555555555550000B 

DATA 6612131415161 7200000B , 2 1 225555555555550000B 

DATA 75552324252627300000B.31325555555555550000B 
DATA 3334353637404 1 420000B , 43445555555555550000B 

END 

0800 

0644,1356,1 



APPENDIX II. PROGRAM FOR FREQUENCY DETERMINATION BY FFT 


//BRIAN JOB 

/•ID PS.---5184 ,NAME=BRIAN GILCHRIST 
/•ID CODE=MARCI 

/•ID PLOT= YES , I0REQ=5000 , LINES=5000 , TIME= (3.00). REGI0N=230K 
/•SETUP UNIT=TAI E , R=A3 1 7 , IDs ( WI534F , NORING ) 

// EXEC FORTLDPT , PARM . G0= ’ TIME= 10.00, PLs60 ' , REGION , GO= 280K 
//F0RT.SY3IN DD • 

C IWK MUST BE FOUND IN CORE ON A BOUNDRY DIVISIBLE BY 8 

C INTEGER*# IWK (18050) WILL ONLY WORK HALF THE TIME 

C THIS IS BECAUSE IT IS USED AS A STORAGE AREA FOR A REAL »4 VARIABLE IN FFT 
REAL*8 IWK(9025) 

INTEGERS ARRAY(32128) 

REAL»8 DATA (6000) 

COMPLEX* 16 GAMN 

REAL*4 DMAG(3000),FREQ(200) 

EQUIVALENCE ( DATA , DMAG ) 

NAM2LIST/PRAMS/NR30 , NSNGO , NELT , NELN 
ASSIGN 999 TO KILL 
CALL TPOPIZ( 12) 

C NRGO IS FIRST RECORD NUMBER 

C NSNGO IS FIRST SCAN NUMBER IN RECORD 
C NELT IS ELEMENT NUMBER IN SCAN 
C NELN IS NUMBER OF ELEMENTS TO BE USTD. 

C WARNING - USER MUST COMPUTEALL OF THESE VALUES CAREFULLY 
C PLOTTING SEQUENCE •••***»»»**»*«•••*»* 

NCRV=35 

C GENERATE FREQ. ARRAY (MINUS TO MOVE DOWN THE PAGE) 

DIMENSION SCA(2) ,SCF(2) 

DO 8 NF= 1,200 
8 FR£Q(NF)=1-NF 

C FIRST DRAW FREQ. CCALE 

C START AT ZERO, THEN GO AT 20. PER INCH(FOR 8 INCH AXIS) 

3CF( 1 )=0.0 
SCF(2)=20. 

CALL PLOT(.5,.5,-3) 

CALL CCP5AX(0. ,9. ,'LOW FREQ(HZ) ' ,-13 ,2. ,-90. ,SCF) 

3CF( 1 ) =380 

CALL CCP5AX(0. ,7. , 'FREQUENCY (HZ) ' ,-1 4 ,7 . , -90 . , SCF) 

C READY SCALE FACTOR FOR FREQ ARRAY 

SCF( 1 ) =— 1 30 . 0 

C SET UP TIME AXIS. (NCRV = # OF PLOTS TO MAKE) 

SCA( 1 ) = 1. 

SCA(2)=1 . 

DIST = NCRV 

CALL CCP5AX(0. ,0. , ' TIME (RELATIVE) ' , -14,DIST,0. ,SCA) 

CALL CCP5AX(0 . , 9 . , ' TIME ( RELATIVE) ' , 14, DIST, 0. ,SCA) 

C PLOTTING SEQUENCE *»»»*•»**•»»**»•*»*• 

4 3 CONTINUE 

READ (5, PRAMS) 

WRITE (6,889) NRGO , NSNGO , NELT , NELN 



889 F0RMAT(1X,9I10) 

C FORMAT DATA CARD &PRAMS WRGO=400O,NSNGO=2O, ETC ,&END 

IF(NRGO.EQ.O) GO TO 42 
C TAPE HANDLING ROUTINES 

IF(NRGO.EQ.I) CALL TPBSRZ( 12) 

IF ( NRGO . EQ . 1 ) NSNGOsKNSNGO 
IF(NRGO.EQ.I) GO TO 41 
KSKIP = NRGO -1 
DO '40 N= 1 , KSKIP 
CALL TPFSRZ( 12) 

40 CONTINUE 

41 FNTOP=FLOAT ( NELN+NSNC?0 ) 

NTOPsIFIX( FNT0P/400 . )+1 
KNSNGO= IFIX ( FNTOP ) +400» ( 1 -NTOP ) 

WRITE ( 6 , 99 ) FNTOP ,NRGO , NSNGO , NTOP , KNSNGO 
99 FORMAT ( * DEBUG * ,F15.6,4I10) 

DO 50 N=1 t NTOP 

M= (N-1 )*2008 + 1 

CALL TPGETZ (12, ARRAY ( M) ) 

CALL TPCHKZ(12,NTAPE,KILL) 

C IGO=M+1999 

C ISTOP = IGO+ 8 

C WRITE(6,104)(ARRAY(I),I=IGO,ISTOP) 

104 F0RMAT(1X,9I10) 

C NTAPE = # OF BYTES TRANSFERRED 

50 CONTINUE 

WRITE(6 ,102) ARRAY(2006),ARRAY(2007),ARRAY(2008) , NTAPE 
102 FORMAT (1 HO,’ TIME CODE ’,4110) 

42 N005=2005 
NA=NSNG0*5+NELT 
DO 10 N=1,NELN 
DATA(N)= ARRAY(NA) 

NA=NA+5 

IFCNA.LE. N005 ) GO TO 10 
NA=NA+8 

N005 = N005+2008 
10 CONTINUE 

TOTAL=0 .0 
DO 26 N=1,NELN 

26 TOTAL=TOTAL+DATA ( N ) 

TOTAL=TOTAL/NELN 

DO 27 N= 1 ,NELN 

27 DATA ( N ) rDATA ( N ) -TOTAL 

25 CONTINUE 

C DC NOW REMOVED 

C AT THIS POINT, "DATA" IS READY FOR WINDOWING 

NELN2=NELN/2 
FLN2=NELN2 

STT=1.f1. /FLOAT (NELN) 

GBASE = 0.01111 
DO 20 K= 1 ,NSLN2 
DEL = 3 . 0* ( STT-FLOAT ( K ) /FLN2 ) 

G=EXP ( -0 . 5*DEL*DEL ) -GBASE 
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DATA ( K ) = G*DATA(K) 

KTP=NELN-tC+1 
DATA(KTP) = G*DATA(KTP) 

20 CONTINUE 

C ######*»«# 

CALL FFTR ( DATA , GAMN , NELN , IWK) 

AT THIS POINT, "DATA " CONTAINS THE COMPLEX SPECTRUM. 
WRITE IT OUT 
CONVERT TO MAGNUTUDES. 

DO 30 N=1,NELN2 
NT = N*2 

DMAG ( N ) =DSQRT ( DATA ( NT ) *DATA ( NT ) +DATA ( NT- 1 ) *DATA ( NT- 1 ) ) 

30 CONTINUE 

C WRITE IT OUT IN SELECTED GROUPS »•*»•»*•* 

DO 31 N=1, 70,10 
NFREQ = N-1 
NTP = N + 9 

31 WRITEC6 , 103) NFREQ, ( DMAG ( K) ,K=N, NTP) 

103 FORMAT ( 1X,I9,10(1X,E9.2)) 

DO 32 N=401 ,551,10 
NFREQ= N-1 
NTP = N + 9 

32 WRITL(6 , 103) NFREQ, ( DMAG ( K) ,K=N, NTP) 

C #»*»«* FIN D FREQUENCY OF PEAK****** 

K 1 = 3 
K2=16 
MAX=5 
53 AMAX=0 . 

DO 52 K=K1 ,K2 

IF(DMAG(K) .LT.AMAX) GO TO 52 
MAX=K 

AMAX=DMAG(K) 

52 CONTINUE 

C COMPUTE FREQ. OF SIG. ALOGORITHIM 
MAXM=MAX-1 
MAXP=MAX+1 

RR = ALOG ( AMAX/DMAG ( MAXP ) ) / ALOG ( AMAX/DMAG ( MAXM) ) 

FMAX1 =MAX-1 

FSIG=FMAX 1 +0 . 5* ( 1 • 0-RR ) / ( 1 . 0+RR ) 

WRITE (6,1 05 )FSIG 

106 FORMAT (/,' FREQ. OF PEAK= * , F 1 0 . 5 ) 

K1 =5 
K2=70 

XKM1 =DMAG( 1 ) 

XKM=DMAG(2) 

XK=DMAG(3) 

XKP=DMAG(4) 

GO TO 6 
5 K 1=400 
K2=540 

XKM1 =DMAG( 39'J ) 

XKM=DMAG(397) 

XK=DMAG( 39 5 ) 
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XKP=DMAG( 399) 

6 DO 54 K=K1,K2 

IF(XK.GT.XKM. AND.XKM.GT.XKM1 .AND.XK.GT.XKP) GO TO 55 
56 XKM1=XKM 
XKM=XK 
XK=XKP 
XKP=DMAG(K) 

GO TO 54 
55 MAXM=K~3 

RR=AL0G ( XK/XKP ) /ALOG (XK/XKM) 

FSIG=MAXM+0 . 5* ( 1 . -RR ) / ( 1 . +RR ) 

WRITE (6,106) FSIG 
GO TO 56 
54 CONTINUE 

IF(K1 .LT.200) GO TC j 
DANGER<<<<<DATA CAN (]) GIVE DIVIDE BY 0 
END OF FREQ. COMPUTATION ON SPECTRAL PEAKS 
PLOTTING SEQUENCE »»»****#»»*#»*«»#»*# 

SET UP SCALING FOR ALL PLOTS 

WARNING*****XLOGZ EXPECTS REF. PT. IN LOWER LEFT<<<<<<<<<<<<<<< 

WARNING DMAG SCALE HAS BEEN SET TO "NICE" VALUES< <<<<<<<<<<< 
WARNING MINIMUM ALLOWED VALUE OF DIST" IS 2. 0<<<<<<<<<<<<<< 
SCA WILL NOW BE SET TO FIXED VALUES... 

SCA ( 1 ) ™ 1 « 

IF(NELT.EQ.5) SCA(1)=0.0 

SCA(2)r.2. 

CALL XL0GZ( DMAG , FREQ , 40 , 1 , SCA , SCF ) 

WRITE(6,101) SCAO) ,SCA(2) 

101 FORMAT ( 1 X , 2E 1 5 . 5 ) 

CALL XLOGZ( DMAG (381 ) ,FREQ(41 ) , 139, 1 , SCA, SCF) 

CALL PLOT(1 . .0. ,-3) 

READY FOR NEXT XLCGZ CALL. 

PLOTTING °EQUENCE 

READ NEXT DATA CARD IF SCAN # = - Z , DO OTHER NELT 

...IF CCAN# = 0, READ NEXT NELN ELEMENTS. 

GO TO 43 
999 WRITE (6,998) NTAPE 
998 FORMAT ( 1 TAPE ERROR ' ,110) 

STOP 

END 

//GO.FT12F001 DD UNIT=TAPE,VOL=( , RETAIN, SER=WI534F) ,DISP:OLD, 

// DCB=(RECFM=U,BLKSIZE=40 16), LABEL: (1,BLP, ,IN) 

//GO.SYSIN DD » 

&PHAMS NRGO= 1 287, NSNGO= 151, NELT= 1 ,NELN=5000,4END 
&PRAMS NRGO= 1 ,NSNG0=1 ,NELT=1 ,NELN=5000 ,&END 
&PRAMS NRGOsI ,NSNG0=1 ,NELTs1 ,NELN-=5000 ,&END 
&PRAMS NRG0=1 ,NSNG0=1 ,NPLT=1 ,NELN=5000,&END 



APPENDIX III. PROGRAM TO CALCULATE ELECTRON CONCENTRATION 
FROM FARADAY ROTATION 


//BRIAN JOB 

/•ID PS=5 1 84 , NAM2=BRI AN GILCHRIST 

/•ID CODE=MARCI 

/•ID TIME=2 

/•ID LINES=4500 

// EXEC FORTLDGO 

//F0RT.3YSIN DD • 

C ANALYSIS OF FR WITHOUT DA, GIVEN CF MODEL. 
C ENTER PARAMETERS OF OBLATE EARTH 
A=5.378l65E6 
A2 =A*A 


C 

C 

c 

c 

c 

c 

c 

50 


C 

c 


A4 = A2*A2 
FLATi1.-1./298.3 
B2=(A*FLAT)**2 
A2B2=A2*(1.-FLAT«*2) 

A4B4=A4*( 1 .-FLAT**4) 

WRITE( 11,41) 

41 FORMAT! '-' , 'NOTE THE FOLLOWING ABREVIATIONS AND CONVENTIONS: V 
*2X , 

*' AZD - LOOK AZIMUTH ------ 

*'ELD - LOOK ELEVATION ------ 

*’HT - HEIGHT - - 

*'V - VELOCITY - -- -- -- - 

• ' RLATD - ROCKET LATITUDE- ------- DEGREES'/ 

• ’ RLNGD - ROCKET LONGITUDE ------- DEGREES'/ 

••FRE - FARADAY ROTATION , EXPERIMENTAL - DEGREES/SECOND'/' 

•'DAE - DIFFERENTIAL ABSORBTION, EXPERIMENTAL - DECIBLES’/’ 

*'P - PRESSURE PASCALS’/' ') 

ENTER PARAMETERS 


• - DEGREES'/’ 

’ ,2 X 

• - DEGREES'/' 

' ,2X 

■ - - METERS'/’ 

’ ,2X 

METERS/SECOND’/’ 

t 

t 


,2X 

,2X 


RKTNO - ROCKET NUMBER 

TM - TIME 

F - FREQUENCY (HERTZ) 

CFM - COLLISION FREQ. MODEL PARAMETER (HERTZ/PASCAL) 

READ (5,50) RKTNO, TM,F, CFM 
FORMAT (F10.4,F10.1,F1C.O,F10.0) 

WPITE(11 ,42) 

42 FORMAT INPUT DATA DECK:') 

WRITE! 11,43) RKTNO, TM , F , CFM 

43 FORMAT ( * - ' ,F7 .4 ,2X ,F6. 1 ,2X, 1PE10 .3 ,2X , 1PE10.3) 

WRITE (11, 40) 

40 FORMAT ( ' ’ ,2X, 'AZIMUTH ELEVATION ALTITUDE VEL. LATITUDE', 
*' LONGITUDE FRE. DAE.') 

CALL COEFF (TM) 

w = 6.2831853*? 

DTR = 1.745329E-2 
INITIALIZE ELECTRON DENSIP 
SD = 1.E9 
ENTER VARIABLES 
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C AZD - AZIMUTH (DEGREES) 

C ELD - ELEVATION (DEGREES) 

C HT - HEIGHT (METERS) 

C V - VELOCITY (METERS/SECOND) 

C RLATD - ROCKET LATITUDE (DEGREES) 

C RLNGD - ROCKET LONGITUDE (DEGREES) 

C FRE - FARADAY ROTATION, EXPERIMENTAL (DEGREES/SECOND) 

C DAE - DIFFERENTIAL ABSOBPTION, EXPERIMENTAL (DECIBELS) 

P - PRESSURE (PASCALS) 

1 RE AD ( 5 , 5 1 , END= 30 ) AZD , ELD , HT , V , RLATD , RLNGD , FRE , DAS 
51 FORMAT (8F10.4) 

WRITE ( 1 1 , 44 ) AZD , ELD , HT , V , RLATD , RLNGD , FRE , DAE 

44 FORMAT ( ' ' ,F9.4,2X,F9.4,2X,F7.0,2X,F5.0,2X,F8.4,2X,F9.4,2X,F8.3,2X 
»,F7.4) 

DAE=0 .0 
READ (5 ,52)P 
52 FORMAT (F10.7) 

WRITE( 11,45) P 

45 FORMAT ( • ’ ,2X,1PE9.2) 

CF = CFM*P 

C CONVERT DEGREES TO RADIANS 
AZ = AZD*DTR 
EL = ELD*DTR 
RLATR = RLATD*DTR 
RLNGR = RLNGD*DTR 
SINLA = SIN(RLATR) 

20 SINLA2=SINLA*SINLA 
C0SLA2=1 .-SINLA2 
SPH = SIN(RLNGR) 

CPH = COS (RLNGR) 

C FIND GEOCENTRIC COORDINATES OF ROCKET 
DEN2=A2-A2B2«3INLA2 
DEN=SQRT(DEN2) 

FAC= ( ( ( HT*DEN ) +A2 ) / ( ( HT»DEN ) +B2 ) ) »*2 
CT=SINLA/SQRT(FAC*COSLA2+SINLA2) 

F=SQRT(HT*(HT+2 .*DEN)+( A4-A4B4*SINLA2) / DEN2) 

21 ST=SQRT(1 .-CT»«2) 

C CALCULATE GEOMAGNETIC FIELD AT ROCKET 

CALL FIELD ( R ,ST ,CT ,SPH ,CPH , BR , BT , BP ,B) 

C TRANSFORM FIELD COMPONENTS, GEOCENTRIC TO GEODETIC. 
SIND=5INLA*ST-SQRT(C0SLA2) # CT 
C0SD=SQRT(1.-SIND*»2) 

BN=-BT»COSD-BR»SIND 

BD=BT*SIND-BR*COSD 

BW=-BP 

S = -1 .758803E1 1 *B 
Y = S/W 

C CALCULATE COSINE OF PROPAGATION ANGLE 
CEL = COS (EL) 

CTH = (CEL»COS(AZ)»BN-SIN(£L)*BD-CEL*S: !(AZ)»BW)/B 
TH=ARCOS(CTH)/DTR 

C CALCULATE FR AND DA COEFFICIENTS 
FV = F*V 
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FC = 6.004153E-7 # FV 
AC = 1 .820428E-7»FV 
C ITERATE TO MATCH FR 
WRITE (6,67) 

67 FORMAT (1H1) 

WRITE (6,68) RKTNO , TM , CFM 

68 FORMAT (3X ,F10 . 4 ,F10 . 1 , 1PE10 .2 , /) 

1 = 6 

DO 2 N=1 ,1 

C CALCULATE REFRACTION AND ABSORPTION INDICES 

CALL SENWYL (ED,CF,3,CTH,W,RI0,RIX, AIO,AIX,ROR,RXR,ROI,RXI) 
FR = FC*(RIO-RIX) 

DA=AC*ABS(AIX-AIO) 

X = 3182.6018«ED/W/W 
Z=CF/W 

WRITE (6, 61 )AZD,BN,ROR,RIO,F,HT 
WRITE (6,62) ELD , BD , RXR , RIX , FRE , ED 
WRITE ( 6 , 63 ) RLATD , BW , ROI , AIO , FR , X 
WRITE ( 6 , 64 ) RLNGD , B , RXI , AIX , DAE , CF 
WRITE(b.65)V,Y,CTH,TH,DA,Z 


61 

FORMAK ' 

AZ 

* , 1 PE 13.6, r 

BN 

' .1PE13.6, ' 

ROR 

' , 1 PEI 3-6, 


& 

RIO 

' , 1PE13-6, ' 

F 

* , 1PE13-6 , ’ 

HT 

* , 1 PE 13.6) 

62 

FORMAT (' 

EL 

* , 1PE13.6, • 

BD 

• , 1PE13-6 , ' 

rx.; 

• , 1PE13-6 , 


& ’ 

RIX 

< , 1PE13-6, • 

^RE 

' , 1 PE 13-6, * 

ED 

• .1PE13.6) 

63 

FORMAT (• 

LAT 

• , 1PE13.6,’ 

BW 

• .1PE13.6, • 

ROI 

' . 1FE1 3 . 6 , 


& 

AIO 

' , 1PE1 3.6,* 

FR 

' .1PE13.6,' 

X 

' , 1PE1 3-6) 

64 

FORMAT ( * 

LNG 

' , 1PE13.6 , » 

B 

* , 1PE13-6, • 

RXI 

' , 1PE1 3 -6 , 


& 

AIX 

• , 1FE1 3 . 6 , * 

DAE 

• ,1PE13-6, ' 

CF 

• ,1PE13.6) 

65 

FORMAT ( ' 

V 

' , 1PE13-6, ' 

Y 

' , 1PE13-6, • 

COS 

* , 1 PE 13-6, 


& 

TH 

' , 1 PE 13.6, ' 

DA 

' , 1PEI3-6, ' 

Z 

' , 1 PE 13-6) 


WRITE (6,60) 

60 FORMAT (1H ) 

ED = FRE/FR*ED 
2 CONTINUE 

GO TO 1 

30 WRITE( 11,40) 

STOP 

END 

SUBROUTINE COEFF(TM) 

COMMON /C0EFFS/TG(11,11) 

DIMENSION G( 1 1 , 1 1 ) ,GT( 11,11) ,SHMIT( 11,11) 

C READ SPHERICAL HARMONIC COEFFICIENTS 

C THE COEFFICIENTS IN THE DATA STATEMENT ARE GSFC( 12/66) 

C THIS G ARRAY CONTAINS BOTH G AND H VALUES 

DATA G/ 

ZO. ,-30401 .2,-1540. 1 ,1307.1,949.3,-233.5,49.2,72.2,8.5, 10. 4,-2. 9, 
15778. 2, -2163. 8, 2997. 9, -1988. 9, 803. 5, 355. 7, 57. 5, -53. 7, 6. 5, 5. 8, -.9, 

2- 1932., 202. 9, 1590. 3, 1276. 8, 502. 9, 228. 4, -.8, 7. 9, -9. 3, 7. 5, -2. 2, 

3- 425.4,227.8,-133.8,881.2,-397.7,-28.8,-238.3,15.6,-9.6,-15.1,. 3, 
4160.3,-274.3,2.3,-246.6,266.5,-157.9,-1.5,-24.3,-6.1,12.1,-2.8, 
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55.1 ,117.8,-114.8,-108.9,82.4,-62.2,-2. ,-3-6, 5. 5, *». 7, 6. 4, 

6- 12. 1,104. 4, 56. 6, -23. 4 ,-14. 8, -13. 3, -108. 9, 15. 5, -8.1,. 2, 4. 7, 

7- 53- 7, -27. 4, -8. 1,7., 24. 3, -22. 5, -21. 4, 3. 6, 13., 1.6, -.2, 

85. 4, -11. 7, 4. 2, -15. 3, 4. 6, 21. 9, -.7, -17. 1,7. 4,. 9, 1.8, 

9-22. 4, 13. 8 , 6 . 3, -3., -1.9, 9., 11. 5,. 1,-1 .5,. 2, 2., 

D-. 1,4. 5,-1., 2. 6, -4. 4,-1 .3, -3-6, 4. , 1 . ,-2. , 1 . 1/ 

DATA GT/ 

Z0., 14. 03, -23-29, -.93, 1-45, 1.61, -.42, -.57,. 35, -.10, -.01, 

1- 3-71, 8. 76, -.09, -10. 62,. 9,. 6,. 82, -.34,. 5, -.13, -.13, 

2- 14. 31 ,-16.62, -4. 56, 2. 31, -1.75, 3. 34,. 82, -1.44, 1.7, -1.2,. 88, 
35. 2, 2. 53, -6. 98, -5. 89,. 66, -.04, 2. 35, -.9, -.11,. 08,-. 18, 

4-2. 19, -.14, 1.88, -6. 52, -3. 01, -.6,. 83,. 03,. 34, -.08,. 17, 

52. 24, 1.59, -2. 61,. 5, -.12, 1.76,. 01, -.6, -.07, -.39, -.02, 

6. 05.. 09. 2. 55, -1.19,. 33,. 84,. 23, -.17,. 43, -.36,. 05, 

7- . 96,. 01,. 43,. 75, -.33,- 49,. 9, -.64, -.15,. 47,. 17, 

8- . 5, -. 21,. 03, -.79,. 05,. 1,-. 36, -.43, -.42,. 37,. 16, 

9. 66.. 54.. 03.. 35, -.03, -.01,. 45, -.35,. 75, -.46,. 31, 

D-. 61, -.64,. 02,. 05, -.63, -.07,. 07, -.03, -.02, -.45, -.23/ 

C CALCULATE NORMALIZATION CONVERSION FACTORS 
SHMIT( 1 , 1 )=-1 . 

MAXN =11 
DO 15 N=2,MAXN 

SHMIT ( N , 1 ) =£HM1T ( N- 1 , 1 ) * FLOAT ( 2 *N-3 ) /FLOAT ( N- 1 ) 

SHMITO ,N)=0. 

JJ=2 

DO 15 M=2,N 

SHMIT (N,M)=SHMIT(N,M-1 )*SQRT( FLOA T( CN-M+1 )*JJ) /FLOAT ( N+M-2 ) ) 
SHMIT (M-1 ,N)=SHMIT(N,M) 

15 JJ=1 

C CONVERT COEFFICIENTS, SCHMIDT TO GAUSS 
DO 16 N=2,MAXN 
DO 16 M=1,N 

G(N,M)=G(N,M)*SHMIT(N,M) 

GT(N,M)=GT(N,M)*SHMIT(N,M) 

IF (M.EQ.1) GO TO 16 

G(M-1 , N ) =G ( M- 1 ,N)*SHMIT(M-1 ,N) 

GT(M-1 ,N)=GT(M-1 ,N)*SHMIT(M-1 ,N) 

16 CONTINUE 

C CONVERT COEFFICIENTS TO NEW TIME 

17 T=TM- 1960.0 

DO 18 N= 1 ,MAXN 
DO 18 M=1,N 

TG(N,M)=G(N,M)+T«GT(N,M) 

IF (M.EQ.1) GO TO 18 

TG(M-1 ,N)=G(M-1 ,N)+T*GT(M-1 ,N) 

18 CONTINUE 
RETURN 
END 

SUBROUTINE FIELD (R,ST,CT,SPH,CPH,BR,BT,BP,B) 

DIMENSION PC 1 1 , 1 1 ) , DP (11,11), CONSTC 11,11) 

DIMENSION SP(11), CP ( 1 1 ) , FN ( 1 1 ) , FM ( 1 1 ) 

COMMON /COEFFS/G( 11,11) 

NMAX = 11 

PCI ,1) = 1. 


1 


A 
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DP ( 1 , 1 ) = 0. 

SP(1) = 0. 

CPO) = 1. 

DO 2 N = 2,11 
FN(N) = N 
DO 2 M = 1 ,N 
FM(M) - M-1 

2 CONST(N,M) = FL0AT( (N-2)**2-(M-1 ) **2) /FLOAT ( (2*N-B)*(2*N-5) ) 

3 SP(2) = SPH 
CP(2) = CPH 

DO 4 M = 3.NMAX 

SP(M) = SP(2)*CP(M-1 )+CP(2)*SP(M-1 ) 

4 CP(M) = CP(2)»CP(M-1)-SP(2)»SP(M-1) 

A0R=6 .3712E6/R 

AR = A0R* # 2 
BT = 0. 

BP = 0. 

BR = 0. 

DO 8 N = 2.NMAX 
AR = A0R*AR 
DO 8 M = 1,N 
IF(M.EQ.N) GO TO 5 
P(N,M)=CT*P(N-1 ,M) 

DP(N,M)=CT»DP(N-1 ,M)-ST*P(N-1 ,M) 

IF(M.EQ.N-I) GO TO 7 
P ( N ,M) =P ( N , M) -CONST ( N , M) *P ( N-2 , M) 

DP(N,M)=DP(N,M) -CONST ( N , M) *DP ( N-2 , M) 

GO TO 7 

5 P(N,N) = ST # P(N-1 ,N-1 ) 

DP(N , N) = ST*DP(N-1 ,N-1 )+CT*P(N-1 ,N-1 ) 

7 PAR = P(N,M)*AR 

IF (M.EQ.1) GO TO 9 

TEMP = G(N,M)*CP(M)+G(M-1 ,N)*SP(M) 

BP = BP-(G(N,M)*SP(M)-G(M-1 ,N)*CP(M) )*FM(M)*PAR 
GO TO 10 

9 TEMP = G( N ,M)*CP(M) 

BP = BP-(G(N ,M)*SP(M) )*FM(M)*PAR 

10 BT = BT+TEMP*DP ( N , M) *AR 

8 BR = BR-TEMP*FN(N)*PAR 

BP = BP /ST 

B = SQRT(BT*BT+BP*BP+BR*BR) 

GTT= 1 .E-9 

BR=BR»GTT 

BT=BT*GTT 

BP=BP*GTT 

B=B*GTT 

RETURN 

END 

SUBROUTINE SENWYL( ED , CF , S , CPH , W , RIO , RIX , AIO , ATX , ROR , RXR , ROI , RXI ) 
C GENERALIZED MAGNETO-IONIC THEORY, SEN-WYLLER EQUATIONS. 
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COMPLEX EI,EII,EIII,AA,BB,CC,DD,EE,VC,VD,VE,VF,VG,VH,VO,VX 
COMPLEX CMPLX , CSQRT , BBVA , CCCPH , RO , RX 

C DEFINE BURKE-HARA FUNCTIONS FOR C3/2 AND C5/2. 

C32(X) = (X»(X»(X»(X+2M. 6531 15)+1 13. 9**16)+1 1.287513) + . 023983 i »74)/( 
1X»(X»(X«(X»(X»(X+24.6568l9)+120.49512)+289.58085)+l49.21254)+9.38? 
27372)+. 018064128) 

C52(X) = (X*(X*(X+6 .6945939 ) +16. 90 1002 )+1 . 1630641 )/(X*(X # (X*(X*(X+ 
1 6. 63 1 4497 )+35. 355257 )+68. 920505) +64. 093464) +4. 3605732) 

P2 = 3182.6018*ED 
WPS = W+S 
WMS = W-S 
Q = P2/W/CF 
R = Q/CF 
T = 2.5*Q 
A1 = W/CF 
A2 = WMS/CF 
A3 = WPS/CF 
A = R»W*C32(A1) 

B = T*C52(A1) 

C = R»WMS*C32(A2) 

D = T*C52(A2) 

E = R»WPS«C32(A3) 

+ F = T*C52(A3) 

El = CMPLX (1.-A.-B) 

£11 = CMPLX (.5*(F-D),.5*(C-E)) 

EIII = CMPLX ( A- . 5* ( C+E ) , B- . 5* ( F+D ) ) 

AA = 2.*EI*(EI+EIII) 

BB = EIII*(EI+EIII)+EII**2 

CC = 2.*EI*EII 

DD = 2 .*EI 

EE = 2.*EIII 

VB = CPH*CPH 

VA = 1 .-VB 

VC = BB*BB»VA»VA-CC»CC*VB 
VD=CSQRT(VC) 

BBVA=BB*VA 
VE=AA+BBVA 
VF = DD+EE*VA 
VG = (VE+VD) /VF 
VH = (VE-VD)/VF 
VO=CSQRT(VG) 

VXrCSQRT(VH) 

C SEPARATE INDICES 
RIO ■ REAL (VO) 

AIO = -AIMAG(VO) 

RIX = REAL (VX) 

AIX = -AIMAG(VX) 

C CALCULATE WAVE POLARIZATIONS 
CCCPH=CC*CPH 
R0=-( (BBVA-VD)/CCCPH) 

ROR=REAL(RO) 

ROI=AIMAG(RO) 

RX=-((BBVA+VD) /CCCPH) 

RXR=REAL( RX) 



RXI=AIMAG(RX) 

RETURN 

END 

/ /G0.FT1 1F001 ID SYSOUT =A , DCB= ( RECFM=FA , LRECL= 133) 


//GO.SYSIN DD * 
14.534 1977.7 

2225000.0 

630000.0 



104.8360 

0.190 

80.4478 

89685.0 

1460.0 37.8055 

-75.3202 

13.29 

104.8680 

0.147 

80.4203 

91113-0 

1451.0 37.8049 

-75.3172 

25.79 

104.8776 

0.123 

80.3932 

92532.0 

1442.0 37.8042 

-75.3142 

36.95 
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APPENDIX IV. PROGRAM TO CALIBRATE PROBE CURRENT 


//BRIAN JOB 

/•ID PS=5 1 84 ,NAME=’ GILCHRIST' 

/•ID CODE=MARCI 
// EXEC FORTLDGO 
//FORT.SYSIN DD • 

REAL LNCHSC 

DIMENSION T( 120) ,TIME(600) ,C( 120,5) , CRNT(600) , AVCRNTC 120) , 

1EN( 100) ,UT(600) , RATIO ( 100) ,SEC( 100) 

C RN = ROCKET NUMBER 

C LNCHHR = HOUR AND MINUTE OF LAUNCH (UT) 

C LNCHSC = SECOND OF LAUNCH (UT) 

READ (5, 90) RN, LNCHHR, LNCHSC 
WR ITE ( 6 , 900 ) RN , LNCHHR , LNCHSC 

900 FORMAT ( IX ,F7 ,5X ,I4,5X ,F5 . 1 ) 

C EN = ELECTRON DENSITY 

C READINGS MUST BEGIN AT THE FIRST INTEGER SECOND AFTER CURRENT READINGS START 
READ(5 ,91 ) (EN(I) ,1=1 ,60) 

WRITE (6, 901 ) (EN(I) ,1=1 ,60) 

901 FORMAT (1X,6(2X,E13*4)) 

C T = TIME ON DATA CARDS 

C C = CURRENT ON DATA CARDS ( 2-DIMENSIONAL ARRAY ) 

DO 111 1=1,120 

READ(5 , 100 ,END=150) T(I) , (C(I, J) , J=1 ,5) 

111 CONTINUE 

WRITE(6,902)(T(I),(C(I,J),J=1,5),I=1,120) 

902 FORMAT ( 1 X , F9 • 2 , 5 E 1 5 • 4 ) 

1=121 

C 

C CONSTRUCT SINGLE SUBSCRIPTED ARRAYS OF TIME AND CURRENT 
C 

150 NCARDS=I-1 

C NCARDS = NUMBER OF CARDS (TIME & CURRENT) THAT HAVE 
C BEEN READ AND WILL BE PROCESSED. 

K=0 

C CRNT = CURRENT ON DATA CARDS ( 1 -DIMENSIONAL ARRAY ) 

C TIME = ELAPSED TIME IN 0.1 SECOND INTERVALS 
C UT = UNIVERSAL TIME 
DO 210 1=1, NCARDS 
DO 210 J=1 ,5 
K=K+1 

CRNT(K)=C(I,J) 

TIME(K) = T(1) + (K— 1 ) / 1 0 . 

UT(K) = LNCHSC + TIME(K) 

C THERE SHOULD BE (N-1) OF THE FOLLOWING CARDS FOR A MAXIMUM E.T. OF N MINUTES 
IF(UT(K).GE.60.0)UT(K)=UT(K)-60.0 
210 CONTINUE 
KM = K 
J=0 

MARK = 0 
DO 220 1=1, KM 
C 
C 

c 


FIND INTEGER SECONDS 



o o o o o o 


ITIME=TIME(I) 

TRT=ITIME 

DIF = TIME(I)-TRT 

IF( DIF. GT. 0.01 .AND. DIF. LT. 0.99 )GO TO 220 
J=J+1 

SEC(J) = TIHE(I) 

C AVERAGE CANNOT BE FOUND FOR FIRST AND LAST SECONDS 
IF( J. EQ. 1 )GO TO 219 
IF(I.GT. (KM-9) )GO TO 219 

MULTIPLY CURRENTS BY WEIGHTING FACTORS AND FIND AVERAGE 

TCRNT =0.0 
DO 215 K=1 , 10 
AK = K 

WCRNT = CRNT (1-1 O+K ) # AK 

215 TCRNT = TCRNT + WCRNT 
DO 216 K=1 ,9 
AK = K 

WCRNT = CRNTU+K) * (10.0-AK) 

216 TCRNT = TCRNT + WCRNT 
AVCRNT(J) = TCRNT/ 100.0 
RATIO(J) = EN(J)/AVCRNT(J) 

GO TO 220 

219 AVCRNT(J) = 0.0 
RATIO (J) = 0.0 

220 CONTINUE 
JMX1 = J+1 
SEC(JMXI) = Q.05 

FORM ARRAY FOR PRINT OUT 

WRITE ( 6 , 301) RN , LNCHHR f LNCHSC 
K=0 

XKM = KM 
HAF = XKM/ 2. 

IHAF = HAF 
HF = IHAF 

IF ( ( HAF -HF ) . GT , 0 . 00 1 ) GO TO 225 
N = IHAF 
GO TO 226 

225 N = IHAF+1 

IF((HAF-HF).GT. 0.99)G0 TO 226 
TIME(N*2) = 0.0 
UT(N*2) = 0.0 
CRNT(N*2) = 0.0 
AVCRNT(JMXI) =0,0 
EN(JMXI) = 0.0 
RATIO ( JMX1 ) = 0.0 

226 JJ = 1 
M = 1 
MM = 50 

231 IF( (N*2) .GT. (MM+50) )GO TO 265 
MARK = 10 
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NTWO = N*2 

IF ( NTWO . LE . MM) GO TO 255 
MM50 = MM + 50 

C THIS IS THE PRINTOUT FOR MORE THAN A HALF PAGE BUT LESS THAN A FULL PAGE 
DO 240 II=NTW0,MM50 
TIME (II) = 0.0 
UT(II) = 0.0 
CRNT(II) = 0.0 
240 CONTINUE 

JEND = JMX1 +5 
DO 250 JI=JMX1 , JEND 
SEC(JI) = 0.0 
AVCRNT(JI) =0.0 
EN(JI) = 0.0 
RATIO (JI) =0.0 
250 CONTINUE 
GO TO 265 

C THIS IS THE PRINTOUT FOR LESS THAN A HALF PAGE 

255 DO 260 L=M,NTWO 

IF(SEC( JJ) .EQ.TIME(L) ) GO TO 256 
WRITE(6 ,650) TIME(L) ,UT(L) , CRNT(L) 

GO TO 260 

256 WRITE(6 -651 ) TIME(L) ,UT(L) ,CRNT(L) ,AVCRNT( JJ) ,EN( JJ) , RATIO ( JJ) 

JJ = JJ + 1 

260 CONTINUE 
GO TO 280 

C THIS IS THE PRINTOUT FOR A FULL PAGE 

265 DO 270 L=M,MM 

IF(SEC( JJ) .EQ.TIME(L) )G0 TO 266 

WRITE (6 , 450) TIME(L) ,UT(L) ,CRNT(L) ,TIME(L+50) ,UT(L+50) ,CRNT(L+50) 

GO TO 270 

266 WRITE(6,451) TIME(L) ,UT(L) ,CRNT(L) ,AVCRNT(JJ) ,EN{ JJ) , RATIO ( JJ) , 
1TIMECL+50) ,UT(L+50) ,CRNT(L+50) , AVCRNT( JJ+5) ,EN(JJ+5) ,RATI0(JJ+5) 

JJ = JJ + 1 
270 CONTINUE 

IF(MARK.EQ. 10)GO TO 280 
JJ = JJ + 5 
M = M + 100 
MM = MM + 100 

WRITE(6 ,301) RN , LNCHHR , LNCHSC 
GO TO 231 
280 WRITE(6 ,700) 

STOP 

C 

90 FORMAT ( F7 . 4 , 5X , 14 , F5 . 1 ) 

91 FORMAT(6E13.4) 

100 FORMA T ( F7 • 1 , 5E 1 3 . 4 ) 

301 F0RMAT(1H1 ,10X,' PROBE CURRENT AND WEIGHTED AVERAGES FOR FLIGHT NO. 

1 M,5X, 'LAUNCH TIME = ' , 14 , ' : * ,F4 . 1 , • UT’ ,// ,3X, 'ET' ,4X, 'UT' 

2 ,9X , ' I* , 12X , ' IA* , 1 1X, 'N* , 1 1X. ’N/IA* , 8X , • ET* ,4X , • UT' ,9X , ' I» , 12X , 

3'IA* ,11X,'N' ,11X,*N/IA* ,/) 

450 FORMAT ( 1 X , 2 ( F5 -1,1X),E13-4, 42X , 2(F5 . I, IS), El 3. 4) 

451 F0RMAT(1X,2(F5.1,1X) ,4E13.4,3X,2(F5.1 ,1X) .4E13-4) 

650 FORMAT ( 1X,2(F5.1,1X),E13«4) 
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651 FORMAT ( 1 X, 2 (F5.1 ,1X),4E13.4) 
700 FORMAT ( 1 H 1 ) 

C 

END 

/ /G0.3YSIN DO • 

14.534 0501 0.0 


O.OOOOEOO 

3.6736E08 

7.1275E08 

1.0231E09 

1.7713E09 

1.1056E09 

9.7248E08 

1.2378E09 

1.7004E09 

2.0154E09 

1.4876E09 

1.7399E09 

7.7521E08 

1.3261E09 

7.8465E08 

5.6499E08 

O.OOOOEOO 

8.8611E08 

7.9762E08 

O.OOOOEOO 

9.1819E08 

7.11 16E08 

O.OOOOEOO 

1.9516E09 

1 .0075E09 

1.0114E09 

2.341 1 £09 

1.8292E09 

1.9388E09 

1.9655E09 

1.6734E09 

1.8177E09 

2.1163E09 

1.4362E09 

2.2506E09 

2.1041E09 

1.4892E09 

1.2826E09 

1.6694E09 

2.3264E09 

1.7722E09 

2.1993E09 

2.3582E09 

1.7651E09 

2.3729E09 

3.335TE09 

1.6570E09 

1 .8439S09 

1.8721 £09 

9.8402E08 

2.6688E09 

O.OOOOEOO 

O.OOOOEOO 

1 .8630E09 

2.3049E09 

1.8107609 

2.5580E09 

1.3913E09 

O.OOOOEOO 

1.1381E09 


68.0 

.1156E-07 

. 1 1 15E-07 

. 1073E-07 

. 1050E-07 

.1047E-07 

68.5 

. 1 146E-07 

. 1350E-07 

.15536-07 

.1704E-07 

.2300E-07 

69.0 

. 3160E-07 

.4149E-07 

.5120E-07 

.5896E-07 

.6550E-07 

59.5 

.7339E-07 

.7703E-07 

.7555E-07 

.8069E-07 

.8908E-07 

70.0 

.1031E-06 

. 1086E-06 

.105 IE-06 

•1047E-06 

. 1 122E-06 

70.5 

.1242E-06 

. 1392E-06 

.1602E-06 

.176 IE— 06 

. 1752E-06 

71 .0 

.1765E-06 

. 171 7E-06 

. 1665E-06 

. 1690E-06 

. 1666E-06 

71.5 

. 1587E-06 

.1524E-06 

.1514E-06 

. 1506E-06 

. 1506E-06 

72.0 

.1946E-06 

.2235E-06 

.1590E-06 

.1534E-06 

.1462E-06 
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APPENDIX V. PROGRAM TO TABULATE AND PLOT ELECTRON CONCENTRATION 


PROGRAM GLOBAL ( INPUT , OUTPUT , TRJTRY , TAPE 1 =TR JTRY , PCOUT , TAPE2=PCOUT , 
JECOUT , TAPE3=ECOUT , INTPL , TAPE4=INTPL ) 

REAL T(60) ,ALTI(60) ,TIM(3500) ,PC(3500) ,ALT(3500) ,ECI(2500) 

REAL ZI(2500) 

REAL BB(9),BD(9),CC(9),DD(9) 

C0MM0N/Q/T(60) ,ALTI(60) ,ITCALL 

REAL NO 

REWIND 1 

REWIND2 

REWIND3 

REWIND4 

DATA YES/"Y"/ 

DATA NO/"N"/ 

DATA EEE/"E"/ 

DATA CCC/"C"/ 

WRITE 1000 
READ 1001, A, B 
WRITE 1002 
READ 1003,C,D,CD 
WRITE 1004 
READ 1001, E, EE 
WRITE 1006 
READ 1003,F,G,FG 
WRITE 1008 

C INPUT TRAJECTORY DATA FROM TAPE1 • 

1=0 

10 1 = 1+1 

READ( 1 , 1 009 ) T ( I ) ,ALTI(I) 

IF(EOF(D) 20,10 

C OUTPUT TRAJECTORY DATA ON TERMINAL FOR VERIFICATION. 

20 11=1-1 

WRITE 1010, (TCI) ,ALTI(I) ,1=1 ,11) 

WRITE 1011 
30 READ 1012, H 

IF(H.EQ.YES) GOTO 50 
IF(H.EQ.NO ) GOTO 40 
WRITE 1013 
GOTO 30 
40 WRITE 1014 
STOP 

C INPUT PROBE CURRENT DATA FROM FILE PCOUT. 

50 DO 70 1=1.3500,5 

READf 2 , 1025) TIM(I) , (PC(I+J-1 ) , J=1 ,5) 

IF(E0F(2) ) 80,70 
70 CONTINUE 
WRITE 1026 
STOP 

80 JJ=I-1 
JN=I-5 

C JJ IS THE INDEX OF THE LAST PROBE CURRENT. 

DO 90 1=1. JN. 5 



DO 90 J=1 ,4 

90 TIM(I+J)=TIM(I)+. 1*J 

C ARRAY TT.M IS NOW FILLED WITH THE TIME OF EACH PROBE CURRENT. 
IF ( TIM( 1 ) .LT.T( 1 ) ) GOTO 92 
IF(TIM( JJ ) .LT.T(II) ) GOTO 55 
WRITE 1028,TIM(JJ),T(II) 

STOP 

92 WRITE 1027 ,TIM( 1 ) ,T( 1 ) 

STOP 

55 WRITE 1015, A, B 

60 READ 1012, H 

IF(H.EQ.YES) GOTO 65 
IF(H.EQ.NO) GOTO 120 
WRITE 1013 
GOTO 60 

65 WRITE 1017 
READ 1018, AM 
WRITE 1019 
READ 1020, PH 
WRITE 1021 
READ 1022, TAU 
WRITE 1023 
READ 1020, T1 
WRITE 1024 
READ 1020, T2 
DO 110 1 = 1, JJ 

IF ( (TIM( I) .GE.T1 ) .AND. (TIM( I) .LE .T2) ) GOTO 100 
GOTO 110 

100 IF(PC(I) .EQ.2.0) GfirO 110 

PC(I)=PC(I)«(1 .0+AM«C0S(2.*3.l4l6«(TIM(I)-PH)/TAU); 

110 CONTINUE 

C ASPECT CORRECTION IS NOW COMPLETE, 

ASPC1 = 1 . 0+AM* COS ( 2 . » 3 • 1 4 1 6 » ( T 1 -PH ) /TAU ) 

ASPC2= 1 . 0+AM*COS (2-*3- I4l6*( T2-PH) /TAU) 

WRITE 1029, A, B, AM, PH, TAU, T1 ,T2,T1 ,ASPC1 ,T2,ASPC2 
120 WRITE 1030 

C CONVERT FROM PROBE CURRENT VERSUS TIME TO VERSUS ALTITUDE. 
ALTMAX=0 . 0 
DO 130 1=1, JJ 
CALL TRAJ(TIM(I) ,ALT(I) ) 

IFCALTMAX.LT. ALT(D) MM=I 

IF ( ALTMAX . LT . ALT ( I ) ) ALTMAX= ALT ( I ) 

130 CONTINUE 

WRITE 1031 , ALTMAX, TIM(MM) 

C CONVERSION TO ELECTRON CONCENTRATION FOLLOWS. 

WRITE 1032, A, B 
READ 1033, NUM 
IF(NUM.EQ.I) GOTO 175 
WRITE 1034 
READ 1035, BD(1) 



133 IJK=1 

CC(IJK)=0.0 
WRITE 1036, IJK 
READ 1037 ,BB( 1 ) 

135 1JK=IJK+1 

IF(IJK.EQ.NUM) GOTO 170 
IF(UK.GT.NUM) GOTO 180 
WRITE 1036, IJK, BD(IJK-I) 

READ 1C35,BD(IJK) 

137 WRITE 1039, IJK 
140 READ 1012, H 

IF(H.EQ.EEE) GOTO 160 
IF(H.EQ.CCC) GOTO 150 
WRITE 1040 
GOTO 140 
150 CC(IJK)=0.0 

WRITE 1036, IJK 
READ 1037 ,BB(IJK) 

GOTO 135 

160 WRITE 1041, IJK 
READ 1037 ,BB(IJK) 

WRITE 1042 
READ 1035,CC(IJK) 

WRITE 1043 
READ 1035,DD(IJK) 

GOTO 135 

170 WRITE 1044, IJK , BD ( IJK- 1 ) , ALTMAX 
BD( IJK )=ALTMAX+. 0048 
GOTO 137 
175 CC(1)=0.0 

BD( 1 )=ALTMAX+.0048 
WRITE 1053 
READ 1037 ,BB( 1 ) 

C DISPLAY N/I CALIBRATION FACTORS ON THE SCREEN. 

180 WRITE 1045, A, B 

WRITE 1046, BD( 1 ) ,BB( 1 ) 

IF(NUM.EQ.I) GOTO 195 

DO 190 IJK=2 ,NUM 

IF ( CC ( I JK ) . GT . 1 . 0 ) GOTO 185 

WRITE 1047 ,BD(IJK-1) ,BD(IJK) ,IJK,BB(IJK) 

GOTO 190 

185 WRITE 1048, BD(IJK-1 ) ,BD(IJK) , IJK, BB(IJK) ,CC(IJK) ,DD(IJK) 
190 CONTINUE 

C CONVERSION FROM PROBE CURRENT TO ELECTRON CONCENTRATION. 

195 DO 230 LMN=1 , JJ 
DO 200 N0P=1 ,NUM 
IF(BD(NOP) .GE. ALT(LMN) ) GOTO 210 
200 CONTINUE 
STOP 

210 IF(PC(LMN) . EQ.2.0) GOTO 230 
lr (CC(NOP) .GT. 1 .0) GOTO 220 
PC(LMN)=PC(LMN)*BB(NOP) 

GOTO 230 

220 PC(LMN)=PC(LMN)*BB(NOF)*EXP( (ALT(LMN)-CC(NOP) )/DD(NOP) ) 
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230 CONTINUE 

C PPOBE CURRENTS HAVE BEEN CHANGED TO ELECTRON CONCENTRATIONS. 

C 

C LINEAR ALTITUDE INTERPOLATION FOLLOWS. 

WRITE 1049 

C FIND THE LOWEST-ALTITUDE ELECTRON CONCENTRATION. 

1=1 

240 IF(PC(I) .NE.2.0) G010 250 
1=1+1 
GOTO 240 

250 IZT=ALT(I)*10.+.C91 
RIZT=IZT+1 
J=1 

ZI(J)=RIZT/10. 

C ZI(1) IS THE LOWEST ALTITUDE FOR WHICH AN ELECTRON 
C CONCENTRATION CAN BE COMPUTED. READY TO INTERPOLATE. 

255 IP=I+1 

260 IF(IP.GT.JJ) GOTO 295 

IF( ALT ( IP ) . LT . AL T ( IP- 1 ) ) GOTO 295 
IF ( PC ( IP ) . EQ . 2 . 0 ) GOTO 265 
GOTO 270 

265 IP=IP+1 GOTO 260 
270 IF ( ALT ( IP ; . LT . ZI ( J ) ) GOTO 275 
GOTO 280 

275 I=IP $ GOTO 255 
C COME HERE FOR INTERPOLATION. 

280 AI=ZI( J)-ALT(I) 

BI=PC(IP)-PCC.) 

CI=ALT ( IP ) -ALT ( I ) 

285 EI=BI*( AI/CI) 

ECI'J)=PC(I)+EI 

286 ZI(J+1)=ZI(J)+.1 
J=J+1 

IF( ZI ( J ) . 3T . ALT ( IP ) ) GOTO 290 
AI=AI+. 1 $ GOTO 285 
2°0 I=IP $ GOTO 255 

C COME HERE WHEN INTERPOLATION FOR ASCENDING PORTION OF FLIGHT IS 
C COMPLETED. 

295 IEND=IP-1 
JEND=J-1 

C JEND IS THE INDEX OF THE LAST ECI TO BE PRINTED FOR ASCENT. 

C ECI STANDS FOR INTERPOLATED ELECTRON CONCENTRATION. 

C 

C IEND IS THE INDEX CF THE LAST PC USED IN THE INTERPOLATION ROUTINE 
C FOR ASCENT. IN THIS PORTION OF THE PROGRAM PC IS THE UNINTERPOLATED 
C ELECTRON CONCENTR-' TICK (NOT THE PROBE CURRENT). 

C 

C OUTPUT THE ASCENDING ELECTRON CONCENTRATION TABLE TO TAPE 4. 

WRITE (4 , 1050 ) A ,B ,C ,D,CD , E, EE, F,G , FG 
WRITE(4 , 1051 ) 

WRITE (4, 1052) 

KM--ZI( 1 ) + .001 

KMD=10*KM 

LD=ZI( 1 )*10 .+.001 
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ITEST=10-(LD-KMD) 

C ITEST IS THE NUMBER OF VALUES TO BE PRINTED IN THE FIRST LINE. 

GOTO( 301 , 302 , 303 , 304 , 305 , 306 ,307,308,309,310), ITEST 

301 WRITE(4 , 1061 )KM,ECI( 1 ) $ GOTO 325 

302 WRITE(4,1062)KM,(ECI(J),J=1, ITEST) $ GOTO 325 

303 WRITE(4,1063)KM,(ECI(J),J=1, ITEST) $ GOTO 325 

304 WRTTE(4 , 106 11 )KM, (ECI( J) , J=1 , ITEST) $ GOTO 325 

305 WRITE(4 , 1065) KM, (ECI( J) , J=1 , ITEST) $ GOTO 325 

306 WRITE(4 , 1066)KM, (EC1( J) , J=1 , ITEST) $ GOTO 325 

307 WRITE(4,1067)KM,(ECI(J),J=1, ITEST) $ GOTO 325 

308 WRITE(4 , 1068)KM, (ECI( J) , J=1 , ITEST) $ GOTO 325 

309 WRITE(4,1069)KM,(ECI(J),J=1, ITEST) $ GOTO 325 

310 WRITE(4 , 1070) KM, (£CI( J) , J=1 , ITEST) $ GOTO 325 
C PRINT ALL BUT THE LAST LINE. 

325 INDEX=ITEST+1 
330 KM=KM+1 

IF( (INDEX+9) .GT. JEND) GOTO 340 
INDEX9=INDEX+9 

WF.ITE(4 , 1070) KM, (ECI( J) , J=INDEX,INDEX9) 

IF (( INDEX+9 ).EQ. JEND) GOTO 350 

INDEX=INDEX+10 

GOTO 330 

C PRINT LAST (PARTIAL) LINE OF DATA. 

340 ITE3T=JEND-INDEX+1 

C ITEST IS THE NUMBER OF VALUES TO BE PRINTED IN THE LAST LINE. 

GOTO( 34 1,342,343, 344 , 345 , 346 , 347 , 348 , 349 ) , ITEST 

341 WRITE(4, 1081 )KM,ECI( JEND) $ GOTO 350 

342 WRITE(4 , 1082)KM, (ECI( J) , J=INDEX, JEND) $ GOTO 350 

343 WRITE(4, 1083)KM, (ECI(J) , J=INDEX, JEND) $ GOTO 350 

344 WRITE (4 , 1084) KM, (ECI( J) , J=INDEX, JEND) $ GOTO 350 

345 WRITE(4 , 1035) KM, (ECI( J) , J=INDEX , JEND) $ GOTO 350 

346 WRITE(4 , 1086) KM, (ECI( J) , J=INDEX, JEND) $ GOTO 350 

347 WRITE(4 , 1087) KM ; (ECI( J) , J=INDEX , JEND) $ GOTO 350 

348 WRITE (4 , 1083) KM, (ECI( J) , J=INDEX , JEND) $ GOTO 350 

349 WRITE(4 , 1089)KM, (ECI(J) , J=INDEX , JEND) 

35C WRITE (4, 1080) 

C TABLE OF ASCENDING ELECTRON CONCENTRATION VALUES IS COMPLETE. 
C234567890 1234567890 1234567890 1234567890 1234567890 1234567890 1234567890 12 
WRITE 1090 , ZI( 1 ) ,ZI( JEND) 

C OUTPUT ELECTRON CONCENTRATION TO FILE 3- 
WRITE(3,1092)A,B,C,D,CD,E,EE,F,G,FG 
WRITEC3, 1 09 ? ) ( ZI ( N ) ,ECI(N) ,N= 1 , JEND) 

C CLOSING REMARKS ARE PRINTED ON SCREEN. 

WRITE 1091 

1000 FORMAT ( " "/" "/" "/" GGGG L" ,8X , "0000" , 4X , ''BBBBB" ,5X , "AAAA" .4X 
•, "L",17X"AER0N0MY'V" G G L",7X,"0 OB BA A L" 
»/" G",8X,"L",7X,"0 0 B B A A L" , 15X , "LABORATORY”/ 

*" G GGG L",7X,"0 0 BBBBB AAAAAA L"/ 

*" G G L",7X,"0 0 B B A A L",19X , "ROCKET"/ 

*" G G L" ,7X , "0 OB BA A L"/ 

*" GGGG LLLLLL 0000 BBBBB A A LLLLLL" , 1 3X, 
•"PROGRAM"/" "/" "/"• PROGRAM GLOBAL PERFORMS THE FOLLOWING TASKS:" 

* /" A. THE PROBE CURRENT PROFILE IS ASPECT CORRECTED IF NEEDED." 

• f " B. TRAJECTORY DATA IS USED TO CREATE A FILE OF PROBE CURRE 


0 . 



•NT VERSUS ALTITUDE."/ " C.N/I CALIBRATION FACTORS ARE USED TO 
% CREATE A FILE OF ELECTRON"/ 5X, "CONCENTRATION VERSUS ALTITUDE."/ 
« " D.A LINEAR ALTITUDE INTERPOLATION ROUTINE IS EMPLOYED TO PRO 

•DUCE A TABLE OF"/ 5X, "ELECTRON CONCENTRATION IN 1/10 KILOMETER I 
•NCREMENTS."/" »/■• ENTER ROCKET TYPE AND NUMBER BELOW.") 

1001 FORMAT (2A10) 

1002 FORMAT (" "/"• ENTER LAUNCH DATE AND TIME.") 

1003 F0RMAT(3A10) 

1004 FORMAT (" "/"• ENTER ZENITH ANGLE.") 

1006 FORMAT (" ■/«• ENTER LAUNCH SITE LOCATION.") 

1008 FORMAT(" "/ "• GLOBAL USES TRAJECTORY DATA FROM FILE TRJTRY. TH 

JESE DATA ARE PRINTED BELOW TO"/ " INSURE THEY HAVE BEEN INPUT 
JTO THE SYSTEM CORRECTLY. (Fj. 1 ,4X,F7.3) n ) 

1009 F0RMAT(F7. 1.F10.3) 

1010 F0RMAT(F5.1 ,4X,F7.3) 

1011 FORMAT ( " "/"• IS THE TRAJECTORY DATA CORRECT AS SHOWN?") 

1012 FORMAT ( 1A1 ) 

1013 FORMAT ("• PLEASE RESPOND WITH YES OR NO.”) 

1014 FORMAT(/"» PLEASE RETYPE THE DATA IN FILE TRJTRY USING THE FORMAT 
J(F5.1,4X,F7.3)") 

1015 FORMAT ( " "/"• PART A - ASPECT CORRECTION."/" IS ASPECT CORRECTION 
J NEEDED FOR ",2A10,"?") 

1017 FGRMAT( /"• ASPECT CORRECTION IS ACCOMPLISHED BY MULTIPLYING THE 

J PROBE CURRENT PROFILE BY"/" THE FUNCTION F(T)=1 .0+AM»COS(2.0«PI» 
J(T-PH)/TAU) . "/ " WHAT IS THE AMPLITUDE CORRECTION FACTOR AM? (X. 

JXX)") 

1018 F0RMATOF4.2) 

1019 FORMAT (" "/"* WHAT IS THE PHASE (PH) IN SECONDS? (XXX. X)") 

1020 F0RMATOF5. 1) 

1021 FORMAT ( " "/"« WHAT IS THE PRECESSION PERIOD (TAU) IN SECONDS? ", 
J"(XX.X)") 

1022 F0RMATOF4.1) 

1023 FORMAT ( " "/"* AT WHAT TIME (SECONDS AFTER LAUNCH) SHOULD ASPECT ". 
J"C0RRECTI0N BEGIN? (XXX. X)") 

1024 FORMAT (/"* AT WHAT TIME SHOULD ASPECT CORRECTION STOP? (XXX. X) n ) 

1025 FORMAT ( F7 . 1 , 5E 1 3 • 4 ) 

1026 FORMAT (" THIS PROGRAM IS DESIGNED TO HANDLE ONLY 3500 PROBE ", 
J"CURRENT VALUES.") 

1027 FORMAT(/"« THE PROBE CURRENT PROFILE BEGINS BEFORE THE FIRST TRAJE 
JCTORY DATA POINT. MORE" /"TRAJECTORY DATA IS NEEDED. THE FIRST PROB 
JE CURRENT OCCURRED AT ",F5.1," SECONDS" /"AFTER LAUNCH. THE FIRST T 
JRAJECTORY DATA POINT IS AT ",F5.1," SECONDS AFTER LAUNCH.") 

1028 FORMAT (/"* THE PROBE CURRENT PROFILE EXTENDS PAST THE LAST ENTRY", 
J" IN THE TRAJECTORY"/"DATA , MORE TRAJECTORY DATA IS NEEDED. THE ", 
J"LAST PROBE CURRENT OCCURRED AT " ,F5. 1 /"SECONDS AFTER LAUNCH. ", 
J"THE LAST TRAJECTORY DATA POINT IS AT ",F5.1," SECONDS AFTER"/ 

J "LAUNCH.") 

1029 FORMATC "/"« O.K. ASPECT CORRECTION FOR ",2A10/" F(T)=1.0+", 

JF4.2,"«COS(2.0«PI*(T-",F5.1,")/",F4.1,")"/" FOR 1 BETWEEN ", 
JF5.1," AND ",F5.1," SECONDS AFTER LAUNCH."/" F(" ,F5. 1 , 
JF5.3/" F(",F5.1, ")--", F5. 3) 

1030 FORMAT(/"» PART B - CONVERSION FROM TIME TO ALTITUDE PROFILE") 

1031 FORMATC IS COMPLETE. APOGEE = ",F7-3," KILOMETERS AT ",F7.3, 

J" SECONDS AFTER LAUNCH.") 
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1032 FORMAT(/"» PART C - CONVERSION FROM PROBE CURRENT TO ELECTRON ", 

^'CONCENTRATION PROFILE."/" HOW MANY N/I CALIBRATION FACTORS ARE " j 

J, "NEEDED FOR ",2A10,"?") 1 

1033 FORMAT (111) 

1034 FORMAT ( /"• O.K. CALIBRATION FACTOR 1 GOES FROM GROUND LEVEL TO ", 1 

J"WHAT ALTITUDE? (XXX. XX)") 

1035 FORMAT( 1F6.2) 

1036 FORMAT(/"* WHAT IS THE VALUE OF N/I (",11,")? (X.XXEXX) OBSERVE 

J"FORMAT") I 

1037 FORMAT (1E7 .2) 

1038 FORMAT (/"* CALIBRATION FACTOR ",I1," GOES FROM ",F6.2," TO 'HIS ", 

J"ALTITUDE? (XXX. XX)") 

1039 FORMAT(/"» SHOULD CALIBRATION FACTOR ",I1," BE A CONSTANT OR.", 

J"EXPONENTIAL FACTOR? (C OR E)") 

1040 FORMAT (” PLEASE COOPERATE BY RESPONDING WITH C OR E.") 

1041 FORMAT (/"* N/I(",I1 ,")=BB , EXP((Z-CC)/DD)"/" WHAT IS BB? ", 

J :, (X.XXEXX)") 

1042 FORMAT(/"« WHAT IS CC? (XXX. XX)") 

1043 FORMAT(/"» WHAT IS DD? (XXX. XX)") 

1044 FORMAT(/"» CALIBRATION FACTOR ",I1," GOES FROM ",F6.2," TO ",F7.3, 

J" KM. (APOGEE).") 

1045 FORMAT (/"* THE N/I CALIBRATION FACTORS FOR ",2A10," ARE:") 

1046 FORMAT ( " 00.00 TO ".0PF6.2,' 1 KM. N/I( 1 ) = ", 1PE8.2) 

* 1047 FORMAT (2X.0PF6. 2," TO ",0PF6.2," KM. N/I(" ,11 , ") =" , 1PE8.2) 

1048 FORMAT (2X.0PF6. 2," TO ".0PF6.2," KM. N/I(",I1 ,")=", 1PE8.2,"*EXP( 

J ( Z-" ,0PF6 . 2, ")/" ,0PF6 .2 , ") ") 

1049 FORMAT( /"* PART D - ALTITUDE INTERPOLATION.") 

1050 F0RMAT(//13X,2A10/13X,3A10/13X, "SOLAR ZENITH ANGLE=” ,2A 10/1 3X, 

J3A10) 

1051 FORMAT (/2X," ALTITUDE !".30X, "ELECTRON CONCENTRATION (CM-3) ASCENT 
J DATA",41X,"!") 

1052 FORMAT (4X,"( KM) !" , 1 13X," ! "/I IX, " ! " ,T17 ,"+.0" ,T28 ,"♦. 1 ",T39, 

J"+.2",T50,"+.3",T61,"+.4",T72,"+.5",T83,"+.6",T94,"+.7",T105, ' 

J"+.8",T116,"+.9",T126,"!"/2X,9("-") ," ! " , 1 13( "-") ! ") 

1053 FORMAT (/"• O.X. CALIBRATION FACTOR 1 GOES FROM GROUND LEVEL ", 

J"TO APOGEE."/" WHAT IS THE VALUE OF i,/I(D? (X.XXEXX) ", 

J"OBSERVE FORMAT") 

1061 FORMAT ( 6X , 13 , 2X , " ! " , 99X , IP1E1 1 .2 , 3X , " ! ") 

1062 FORMAT (6X, 13. 2X,"!",88X, 1P2E1 1 .2,3X,"!") 

1063 FORMAT(6X,I3,2X,"!",77X,1P3E11.2,3X,"!") 

1064 F0RMAT(5X, 13 , 2X , " ! " ,66X , *?4E1 1 .2, 3X, " ! ") 

1065 FORMAT (6X , 13 ,2X , " ! " ,55X , 1P5E1 1 .2 , 3X," ! ") 

1 066 FORMAT (6X , 13 , 2X , " ! " , 44X , 1 P6E1 1 . 2 , 3X , " ! " ) 

1067 FORMAT (6X, 13, 2X," ! " , 33X , 1P7E1 1 .2 , 3X , " ! ") 

1068 FORMAT(6X,I3 ,2X," !" ,22X, 1P8E1 1 .2 , 3X," » ") 

1069 FORMAT(6X ,13 ,2X ! ", 1 1X, 1P9E1 1 .2 , 3X, " ! ") 

1070 FORMAT (6X , 13 ,2X, " ! " , 1P19E1 1 .2, 3X , " ! ") 

1060 FORMAT(2X,124("-"1) 

1081 FORMAT(6X ,I3i2X," ! ", 1P1E1 1 .2,9( 1 IX) , 3X," !") 

1082 FORMAT(6X,I3,2X,"!", 1P2E1 1 .2,8(1 IX) ,3X,"!") 

1083 FORMAT(6X,I3,2X,"J", !P3E11.2,7(11X),3X,"!") 

1084 FORMAT (6X , 13 ,2X , " ! " , 1P4E1 1 .2 ,6( 1 IX) , 3X , " ! ") 

1085 FORMAT (6X , I3,2X ," ! " , 1P5E1 1 .2 ,5( 1 IX) ,3X,"1") 

1086 FORMAT(6X, I3,2X,"!",1P6E11. 2,4(1 IX), 3X,"f") 
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1087 FORMAT(6X,I3,2X,"r ,1P7E11.2,3(11X) ,3X,"I") 

1088 FORMAT!6X, 13, 2X, "!", 11 3E1 1 .2,2(1 IX) ,3X,"!") 

1089 F0RMAT(6X, 13 ,2X , " ! " , 1P9E1 1 .2 , 1 ( 1 IX) , 3X , " ! ") 

1090 FORMAT!" MINIMUM ALTITUDE:", F5. 1/" MAXIMUM ALTITUDE:" ,F5 . 1 ) 

1091 FORMAT! " "/"* PROGRAM GLOBAL IS NOW FINISHED. REMEMBER THAT THE", 
J" TWO OUTPUT FILES ARE"/" LOCAL FILES AND MUST BE SAVED OR ", 
J"R£PLACED. OUTPUT FILE INTPL IS THE"/" FILE WHICH CONTAINS ", 
J"THE INTERPOLATED ELECTRON CONCENTRATION IN CHART"/" FORM. ", 
J"OUTPUT FILE ECOUT CONTAINS ALTITUDE/CONCENTRATION PAIRS IN A", 

J" FORMAT"/" 4(0PF10.2, 1PE10.2) SUITABLE FOR PLOTTING.") 

1092 FORMAT ( 2A 10/3A 10/"SOLAR ZENITH ANGLE=" ,2A10/3A10) 

1093 FORMAT! 'U0PF10 .2 . 1PE10 .2) ) 

STOP 

END 

SUBROUTINE TRA J ( TIME , HEIGHT) 

C COMPUTES INTERPOLATED ALTITUDE VALUE FOR ANY TIME AFTER LAUNCH. 

C MODIFIED FOR USE WITH PROGRAM GLOBAL. ALTITUDE AND TIME V * ' UCS 
C AT 10 SECOND INTERVALS MUST BE PROVIDED VIA ARRAYS T AND ALTI 
C FOR ENTIRE PERIOD OF FLIGHT. GLOBAL READS ISESE ARRAYS FROM FILE 
C TRJTRY. 

C INPUT ARGUMENT: TIME = TIME IN SECONDS AFTER LAUNCH. 

C OUTPUT ARGUMENT: HEIGHT = ALTITUDE IN KILOMETERS. 

C 

COMMON/Q/TC60) , ALTI (60) ,ITCALL 
INCR=0 

C IF FIRST CALL TO TRAJ, INITIALIZE VARIABLES. IF NOT SKIP TO 10. 
IF(ITCALL.GT.O) GOTO 10 
1=3 

ITCALL= 1 
5 INCR=1 
IM1=I-1 
IM2=I-2 

10 IF(TIME.LT.T(IM2) )WRITE 801, TIME 
801 FORMAT!" "/"TIME LESS THAN LOWEST TRAJECTORY POINT, TIME=" , F7 . 3) 
C LOCATE TIME VALUES WHICH BRACKET PRESENT INPUT TIME VALUE. 

IF( TIME .LE .T(I ) ) GOTO 15 
1 = 1+1 
GOTO 5 

C IF PRESENT TIME VALUE IS IN SAME INTERVAL AS PREVIOUS ONE, 

C COMPUTE ALTITUDE WITH OLD COEFFICIENTS. IF NOT, COMPUTE NEW 
C COEFFICIENTS FIRST. 

15 IF(INCR.EQ.I) GOTO 25 
20 HEIGHT=A*TIME*TIME+B*TIME+C 
RETURN 

25 BRAC1 = (T! I )-T( IM1 ) )*( ALTI( IM1 )-ALTI( IM2) ) 

BRAC2=(T(IM1 )-T(IM2) )*(ALTI(I)-ALTI(IM1 ) ) 

T0P=BRAC1-BRAC2 

BRAC1 = (T( IM1 )-T( IM2) ) * (T( I )*T( 1)-T( IM1 )®T(IM1 ) ) 

BRAC2= (T( I )-T(IM1 ) )*(T(IM1 )*T(IM1 )-T(IM2) *T( IM2) ) 

BOTTOM=BRAC2-BRAC 1 

A=TOP/BOTTOM 

B= ( ALTI ( IM1 )-ALTI ( IM2) ) -A* (T( IM1 ) *T( IM1 )-T( IM2 ) *T( IM? ) ) 

B=B/(T( IM1 )-T( IM2) ) 

C=ALTI ( IM2) -A*T( IM2)*T( IM2)-B*T( IM2) 

GOTO 20 
END 


